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EDITOR’S NOTE 


The introduction to the first issue of the Irish Undergraduate Mathematical Mag- 
azine, published 15 months ago, expressed the hope that the magazine could 
become a truly national publication. With the present issue containing articles 
from authors in six of the seven universities of the Republic of Ireland, it may 
be argued that that hope has been realized. 

However, with all those who have been involved organizing the magazine 
in the United States next year, there will be no issue of the UMM in 2014 un- 
less there is a change of the guard. Any undergraduate student (or preferably 
students) in Ireland who feel up to the challenge of running the magazine are 
welcome to it. The only job requirement is a strong belief in the idea of commu- 
nicating mathematics in an accessible way, and at least one student involved 
should have a reasonable knowledge of TEX typesetting. New editors will be 
given all of the existing typesetting files, and control of the domain name and 
email address. 


° 


The present magazine would, of course, not have been possible without the 
work of the contributors. Prof. Anthony O’Farrell generously volunteered his 
article to the magazine following his talk on the topic at the 2012 student con- 
ference in UCD; I am very grateful for his contribution. To the diverse group 
of undergraduate students who contributed to this issue: first my apologies, 
that the magazine took so long to go to print, but most of all my thanks for 
your hard work. I hope that you are happy with the finished product! 

Finally, on a more personal note, I would like to thank to Dr. Tom Carroll 
and Niamh Tobin for each providing encouragement in their own way. 
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AN INTRODUCTION TO THE CALCULUS OF 
VARIATIONS 


KTERAN COONEY 
University College Cork 


In this article I hope to present an introduction to the calculus of variations 
that will not only whet the appetite of the reader for the subject, but that will 
also successfully describe an important mathematical tool that can help solve 
a very general range of problems. Most of us will be used to the “calculus 
of functions,” such as taking derivatives or integrals of x? or sin(x), and from 
these gaining a much better understanding of the original functions. For many 
studying mathematics, the day you learn calculus is the first day of the rest of 
your life. In the calculus of variations, instead of studying functions of num- 
bers, we study functions of functions, and find the analogous derivatives of 
these objects. This subject is far more general then its famous predecessor, but 
unfortunately not quite as developed. 


1. FUNCTIONALS 


Functions are very general things, so we need some rules to govern them be- 
fore we start treating them like numbers. To begin with, we will consider 
functions from the real numbers to the real numbers. The set of continuous 
functions defined on a given closed interval [a,b], which we call F, is a vec- 
tor space over the real numbers in the sense that addition and multiplication 
behave how we'd like them to. So we can add and subtract functions to get 
functions and we can multiply them by numbers to get functions. In general, 
we will consider J,,, the set of all functions which map a given closed interval 
[a, b] to R and which can be differentiated n times. Clearly F, is also a vector 
space over the reals. 

Now, we must introduce a notion of distance between our functions. We 
need this idea of distance, because to define differentiation we need to be able 
to take limits. Taking limits is equivalent to “getting close” to something, so 
naturally we need to know what “close” means. We call this notion of distance 
the norm of an element and we denote it by || ||. We want our notion of anorm 
to satisfy 3 conditions: for x and y vectors and w a scaler (real number) we want 


1. |x| =O = x =0. 
2. |lax|| = lel |lxI- 
3. |x + yl < Ixl + My. 


Condition 3 is the familiar triangle inequality. It turns out that on F the most 
intuitive norm is also the most practical. For a given interval |a, b] we define 


FC) = max O91, (1.1) 
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We leave it as an exercise to the reader to show that this definition satisfies the 
above conditions. 

This norm may appear overly simplistic. When I was studying the topic, I 
was more inclined to think that 


b 
fll = i if@lae 


would be a better norm as it seems to capture the “distance” of a function 
better. Indeed, it works perfectly also. However as we are really looking at 
the limiting behaviour of functions as opposed to their global behaviour, the 
above simpler definition suffices. 

So now, using the first norm on the interval [0, 27], we have 


|| A sin(@x)|| = | Al. 


This makes intuitive sense, as the distance of a function is simply the largest 
distance the function is away from the x-axis. We could extend this definition 
to F,, but there would be a problem. If we want a function to “get close” to 
another function, we want it to become very like that function in our limit. In 
our example above, the function A sin(wx) goes to zero as A goes to zero. But 
in J, if we decrease A but increase w, the norm of A sin(wx) will certainly go to 
zero but the derivative Aw cos(wx) could stay constant or even diverge! In F; 
we are studying functions and their derivatives so clearly this idea of closeness 
is wrong. To compensate for this, given a closed interval [a,b] and an element 
F(x) of F, we define 


n 
= (i) 
Ife = 0 max foo]. 

This definition considers the size of the function and its derivatives so that 
when two functions get close to each other they become alike in the sense we 
are interested in. 

Now that we understand our “numbers” correctly, let’s now examine our 
“functions”. A functional is a function which takes a function and assigns to it 
a number. Let’s take F, for a given [a, 5]. Then 


Jf) = f@), J{A(x)] = IFO), 
b 
J{f(x)] = | feats. and JLf@)] = rorre) 


are all functionals. An interesting example of a functional is the total length 
of the graph of a function, which will be derived later in the article. We say a 
functional J on F,, is continuous at a point (remember that f is a function 
now!) if for every positive e there exists a positive 5 such that 


I — yl <8 > |JP]-J iy] <e. 
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We can see why our norm definition was so important. If we had a functional 

that depended on the derivative of its argument but were still using the first 

norm (1.1) then the above definition of continuity would not make much sense. 
We now introduce the variation of a functional. Set 


AJ[h] = Jly + A] — Jly] 


where y and / are both functions and J is our functional. Here h is behaving 
like our increment, or small difference, as it would in regular calculus. Sup- 
pose that we can write 

AJ[h] = ¢[h] + é[h] 


for some e[h] that goes to zero as ||h|| goes to zero. Then we say ¢ is the varia- 
tion of J. This is our “derivative”. According to our definition, it behaves like 
a linear approximation, exactly as in the regular calculus case. In variational 
calculus, we do not “take derivatives” as much as we would in regular calcu- 
lus, as the definition is too broad. Understanding what the variational is and 
that it exists though is key to the subject. 

The most common type of functionals are of the form 


b b 


7S / FO. 9@),9/@)dx = / F(x, y,y")ax (1.2) 


a a 


where F is some real-valued function. These functionals are useful in that 
they can be constructed to depend on the local behaviour of the graph of a 
function at every point. For example, 


b 


Iy|= i Jit y'@dx 


a 


gives us the arc length of the graph of the function y(x). You may ask the 
question why we say that y and y’ are two independent parameters, as the 
derivative of a function is obviously defined by the function. This notation is 
convenient because it indicates explicitly what F depends on and also because 
we are often able to treat y and y’ independently of each other in variation 
problems. In truth, the notation is just a matter of convenience. 


2. THE EULER-LAGRANGE EQUATION 


With all those definitions out of the way, let’s start to apply what we have de- 
veloped! Suppose that we are given a functional of the form (1.2) over F2 for 
some given interval [a,b]. Also, suppose we are restricted to functions with 
known end-points; that is, we know that y(a) = A and y(b) = B for some 
fixed A and B. We are looking for extremals of the functional J. An extremal 
is a function that is locally a maximum or a minimum of the functional J. 
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Clearly we require that the variation of J at y be zero if we want it to be an ex- 
tremal. The proof of this statement is almost identical to the analogous regular 
calculus proof. We consider the effect of a small increment / on our functional: 


AJ = Jly +h] — Jy] 
b b 


= f Foy thy +i dx f Foy yde. 


a a 


Note that h(a) = h(b) = 0 to satisfy our boundary conditions. We now take a 
Taylor series expansion of F about y. There is no difficulty in doing this as F 
is just a real valued function. Then, taking terms up to first order gives 


b b 
AD = fPO.y.9) + OK I+ Brey. x dx = f FCs yy dx 


a a 


where Fy and Fy denote the partial derivatives of F with respect to y and 
y’ respectively. As we have neglected all but the linear terms of our Taylor 
expansion, AJ is in fact our variation. Using the linearity of integration and 
applying our extremal condition g[h] = 0 we get 


b 
g|h] = [@ (x,y, y)ht+ Fy(x,y,y)h’)dx = 0. 


Now we apply integration by parts to the second term in order to remove the 
dependence on h’, remembering that h(a) = h(b) = 0. This gives 


b b 
i Fy (x, y,y)h'dx = Fy (b)h(b) — Fy (ajh(a) -— / nt Fyds 


a 


; d 
= ~ f rR Fyax 


b 
d 
g|h] = i G - + Fy) dx = 0. 


We need this integral to vanish for any significantly small h. This implies that 


and so then 


d 
Fy - ee 0. (1.3) 
This equation is known as the Euler-Lagrange equation. This is the big equation 
in the calculus of variations. 
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As an example, let’s try and find an extremal of the functional 


m/2 


Jb] = i (y? — y’?) dx 
0 


where y(0) = 1 and y(*/2) = 0. We have F(x, y, y’) = y?—y” so Fy = 2y and 
Fy = —2y’. Substituting this result into the Euler-Lagrange equation gives 


d 
wa I) = Seta 0 


which is just our familiar ordinary differential equation for sinusoidal curves. 
So y = Acos(x) + Bsin(x). By substituting in our boundary conditions we 
see that y = cos(x) is our extremal. The questions remains as to whether 
y = cos(x) maximises or minimises the functional. It can be easily shown that 
it maximises it by looking at another function which satisfies the boundary 
conditions, say y(x) = 1— 2x and showing that it has a smaller value. 

After all those proofs and definitions, let me give you two important state- 
ments without proof which really illustrate the importance of the Euler- 
Lagrange equation. If we have a functional of the form 


b 


J [y1, ¥2, +s Yn] = [ Fe. Vis Voie Sis Ins Vee Vocal 


a 


where each y; is a function of x with fixed end points, then at an extremal we 
must have 


d 
Fy, - anf = 0 for alli. 


So finding a curve in R” that maximises or minimises a functional reduces to 
just solving n different Euler-Lagrange equations! Secondly, the Euler-Lagrange 
equation is invariant under a change of variables. So if we are trying to find 
an extremal curve in 3-dimensions, we must solve 3 different Euler-Lagrange 
equations, but we can choose which co-ordinate system to use. We could use 
Cartesians, cylindrical polars, spherical polars, parabolic, hyperbolic, and so 
on. These two properties make this equation exceptionally general. 

It should be clear to the reader that we haven't approached this topic from 
a strictly rigorous point of view, and there are some intricacies that we have 
simply ignored. This said, many satisfactory and rigorous resources on the 
topic may be found (see [1]). 


3. BUBBLES AND MINIMAL SURFACES 


Mathematicians aren’t happy unless there is a problem, so let’s make one. I 
want to find the general shape of a bubble or soap film with a given boundary. 
For instance, what happens if I dip a helix into a soapy mixture and lift it out 
again. What shape will the soap film make? Intuitively, we know that soap 


2013] CALCULUS OF VARIATIONS 9 


film is stretchy; i.e., there is a tension in the soap film that makes it contract 
as much as possible. This is the same thing as saying that a soap film tries 
to minimise its surface area as far as possible. In this way, we call bubbles 
minimal surfaces, in that they minimise surface area. The task remains to define 
the functional J which describes a given surface area. 


Let’s consider a more specific question. Suppose that I have two concentric, 
parallel rings which are a distance d apart and have radii r; and rz as shown 
in the diagram. If soap film was stretched between the two rings, what shape 
would it take? The shape is going to be a surface of revolution because of the 
axial symmetry of the construction. Suppose that we rotate a function y(x) 
about the x axis. We can break up the graph into a sequence of co-axial rings. 
The perimeter of each ring will then be 277y(x). Notice that if we take the limit 
of this sequence we must also account for the local behaviour of the graph y(x). 
What we actually need is the local length of the graph y(x). We presented an 
equation for this earlier, now let's justify it. 


In the below diagram we have an arbitrary function y(x). Here d/ is the 
length of an infinitesimal part of the graph. By Pythagoras’ Theorem we have 
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dl? = dx? + dy?. Factoring out a dx” and taking square roots we see that 


asia (2 a 
= —]} dx 
dx 


giving us the above equation for the length of the curve y, denoted by L[y]: 


b 
Ly|= f i+ yas. 


From this functional, if you wished, you could use the Euler-Lagrange 
equation to show that the straight line is indeed the shortest path between 
two points on the plane! Going back to the surface area problem, the local 
area of a surface of revolution of y is the perimeter of the loop times the local 
length of the graph. If we denote the surface area by S[y] then 


b 


S[y] = for Vil+y?dx. 


a 


Now that we have our functional, we may apply the Euler-Lagrange equa- 
tion. As the S is independent of x, the Euler-Lagrange equation can be sim- 
plified ([1]) to 

F-y'Fy =C. 


So in this instance the Euler-Lagrange equation becomes 


12 
y ay ee eee =C 


Vit y? 


We multiply across by /1 + y”? and factorise, 
PO PP Yaa CVF", 
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re-arrange terms, 


eae eet Os 
IN ag 


and separate variables to obtain 


An integration yields 


2_C¢C2 
re q=cin(2 AYE (eae). 


The solution for this turns out to be a coshine function, defined by 
1 x —Xx 
cosh(x) = ze +e *). 


This curve is know as a catenary. It is a very familiar shape because any time 
you hang a chord or a chain between two points, the resulting shape is a cate- 
nary (guess what, we can prove this using variational techniques). It is also 
used in constructing bridges and structural arcs as it supports weight very ef- 
ficiently. We define the coshine function by If you rotate a catenary about its x 
axis then the resulting shape is a catenoid. It is the surface a soap film takes be- 
tween two circular rings, and is the solution to our problem. To determine the 
values of our constants of integration C and C; we simply insert our boundary 
conditions. 
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¥ 


The problem of finding a minimal surface given its boundary is known as 
Plateau’s problem. Using a knowledge of differential geometry, we can find an 
equation that describes these surfaces (the surfaces turn out to be those with 
a mean curvature of zero). This is a very current field as the use of computer 
graphics has helped stimulate it enormously. In the 1970’s, a Brazilian gradu- 
ate student Celso José da Costa apparently solved the equations for a unique 
minimal surface in his sleep! The surface is known as Costa’s surface and is very 
interesting in that it has a hole in it. One may find many inspiring images of 
these surfaces online, for instance at [4]. 

All this said about bubbles, we never explained why the regular ones are 
spheres! These bubbles are boundary less, so we need “another” boundary 
condition. We know that if there is a volume of air trapped inside the bubble, 
it must stay there. So the total volume enclosed by the bubble is fixed. As we 
know the bubble tries to minimise it’s surface area, we want to find the shape 
that has the least surface area for a given volume. Intuitively, this is the sphere. 
And guess how we would show this is true? (I'll give you a hint: I’m writing 
an article on it.) 

We may also extend (or contract) the idea of minimal surfaces to minimal 
curves or geodesics. These are curves which when restricted to a given surface 
minimise the distance locally. If we attempt to do geometry on surfaces, these 
are analogous to our lines in plane geometry. Again, we may use differential 
geometry and variational techniques to derive the equations that govern them. 


4. THE BRACHISTRONE CURVE 


We now deal with another variational problem, one of the first actually. Imag- 
ine that we are given two points in space and that we want to connect the two 
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points A and B by acurve. If we drop a ball on the curve at A, the ball will roll 
to B. We want to find the curve such that the ball will get from A to B in the 
least amount of time. This curve is known as the brachistrone curve. You may 
be tempted to say the curve will be the straight line, but it is not. We shall see 
why in a moment. The general velocity of a particle, by definition, is given by 
ds dx 
= —_—_—= 1 12 

ca Par 
where the second equality follows as before. We also know from conservation 
of energy ([3]) that E = mgh + mv?/2 is constant. We will set h = 0 at the 
point A. Remembering that the ball will start from rest, set v = ./2gy. So 


dx JV1+y? 
2gy=v=yvl 12 dt = —— 
V2gy=v=yl+y F > Tay dx 


and the total time taken is given by the functional 


T[y] = =| Go 


By applying the Euler-Lagrange equation again and solving for it as we did 
for the catenary, we obtain a family of solutions of the form 


—-,; (9 — cos(@)). 


1 
x(0) = 702° —sin(@)) + C2, y(@) = 3c 


These curves are known as cycloids. They are the result of rolling a circle on 
a straight line and tracing out the path of a point on the circumference. They 
work because there is a steep initial slant to build up velocity and a gentle 
evening out to conserve it. There is another problem very similar to this one, 
the tautochrone problem. In this instance we wish to find a curve such that a 
ball will reach the bottom of it at some after some constant time T regardless 
of where it is dropped on the curve. Interestingly, this curve turns out to be 
the cycloid again! 


5. ODDS AND ENDS AND THE END 


The Calculus of Variations can be used in a wide range of topics, but excels 
particularly in physics (see [3]). So in the following paragraph, let me blurt 
out what else the calculus of variations allows us to do. 
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If we let V be the potential energy of a system and T be the kinetic energy 
of a body at a point in the system, then we define the Lagrangian by L = T—V. 
It turns out (see [2]) that Newton’s laws of physics are equivalent to saying that 
the Lagrangian is minimised. This formulation is known as the principle of least 
action. Fermat's principle states that light will travel in such a way as to minimise 
the time of it’s journey. Einstein’s theory of relativity states that we do not live 
in a Euclidean universe and that if nothing gets in our way we travel along 
the geodesics of the universe. We may use variational techniques in electrody- 
namics to derive the equations of motion. Finally, Feynman took integrals over 
functionals instead of derivatives, and used this work to re-phrase quantum 
mechanics. 

I hope Ihave opened your mind to a genuinely interesting field, and that I 
send you forth unto the world knowing how to use (or at least appreciate) 


d 
Fy ~ = Fy =0. 
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INTEGRAL BINOMIAL COEFFICIENTS 


PROF. ANTHONY O’FARRELL 
National University of Ireland, Maynooth 


1. INTRODUCTION 


The binomial coefficients are defined by 


()-“e 


k k! 


for nonnegative integer k and any a. Usually, a is a real or complex number, 
but the definition makes sense if a belongs to any field of characteristic zero. 
The following is well-known: 


Theorem 1.1. The binomial coefficients (71) are positive integers, for integers n,k 
with0<k <n. 


The usual proof uses the Law of Pascal’s Triangle, and induction. 


The binomial coefficients (;), with rational r, occur in the Maclaurin series 


expansion of (1 + x)” (convergent for real or complex x with |x| < 1). For 
instance, 


© fi 
Vl+x= a tee 
k=0 
Calculating a few terms, one finds that the series begins 
3. 9 4 


(ck ee 
XxX — XxX — eee 
2* 8 ie” 64. 


The coefficients are not integral (or nonnegative), but when common factors 
are cancelled (i.e. they are expressed in reduced form m/n, withm € Z,n € 
N, and gcd(m,n) = 1), it is remarkable that only powers of 2 occur in the 
denominators. This is not an accident: the pattern continues forever. We have 
the following, slightly less well-known result: 


Theorem 1.2. Letr ¢ Qand0 <k € Z. Suppose that r = m/n in reduced form. 
Then the binomial coefficient (;) has reduced form s/t, where t is a product of powers 
of primes that divide n. 


For instance, in the expansion of (1 + x) é , the coefficients all take the form 
s/(2%3°), for some s € Z. 


The theorem may be proved using elementary number theory, for instance 
by reducing it to the statement that if d,k € N and r is the largest factor of k! 
prime to d, then r divides the product of the terms of each k-term arithmetic 
progression of integers having step d. 
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The purpose of this note is to give a very short soft proof of Theorem 1.2, 
by using a little analysis. Specifically, the proof uses the field Q, of p-adic 
numbers. For the benefit of readers who have not met these numbers, we give 
a short introduction in the next section, and then give the proof in the final 
section. 


2. THE P-ADIC NUMBERS 


For a prime p, the p-adic valuation of a rational number is defined by setting 
|Ol|p = O and 


trp” 
Ss 


P 


whenever r € N,s € N,andn € Zwith gcd(r, p) = gcd(s, p) = 1. For instance, 


1 
300 


1 
[200] = 7, [301 = 1, and | 


2 

Thus some numbers that have large absolute value have small valuations, and 
vice-versa. Also, numbers that have small valuations with respect to one prime 
may have large valuations with respect to another. 

The p-adic metric on the set Q is defined by setting the distance between 
two rationals a and b equal to ||a — b||p. You can verify easily that this does, 
indeed, define a metric. In particular, the triangle inequality follows from a 
stronger form known as the ultrametric inequality: 


lla — bllp S maxilla — ellp, lle — bly}. 


The space Q, of p-adic numbers is the completion of Q with respect to the 
p-adic metric. It is a complete metric field, i.e. the field operations are contin- 
uous. One can show (although we do not need this for the proof below) that 
Q, has the same cardinality as R, and that it is locally-compact and totally- 
disconnected. 

The closure of Z in Q, is denoted Z,, and called the set of p-adic integers. It 
is a compact, totally-disconnected metric space, and an integral domain, and 
Q, is its quotient field. 

From the point of view of number theorists, there is little to choose between 
R and any of the Q,. They are all more-or-less equally-interesting ways to 
complete the set of rationals. For instance, if one is interested in solving a 
Diophantine equation such as x? + y? = z? for integers, then it is necessary 
that the equation have a solution in each Zp and in R. For some equations, the 
converse holds — such a result is called a “Hasse Principle". 

Each infinite series of the form 


oo 
) n 

an P 
n=0 
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with a, € Z is convergent in p-adic metric, and so represents some p-adic 
integer. For instance, in 2-adic metric we have the formula 


14+2+44---+2%+4---=-1, 


which may be found in Euler’s work. More generally, for any prime p, 


(p-1)+(p—- Ip t+ (p- 1)? p++ = 1 


in p-adic metric. From this we deduce that every p-adic integer is the limit of 
a sequence of positive integers. 
A non-integral rational number may be a p-adic integer. For instance, 


1 
bp AS pa 


in 3-adic metric. More generally, it is not hard to see that a rational number r 
with reduced form m/n belongs to Z, if and only if p does not divide n. 


3. THE PROOF 


Theorem 3.1. If p € Nis prime, a € Zp and 0 > k € Z, then ({) € Zp. 


Proof. Fix k € Z,k > 0. The function 


ro) 


is a polynomial with coefficients in Q, and hence it is continuous, as a func- 
tion from Q, into Q,. (This just depends on the fact that Q, is a metric field.) 
Choose a sequence (d;)°2., C N with a, — a in p-adic metric. Then f(an) € 
N Cc Zy, and hence f(a) = limn—+oo f(an) € Zp, since Zp is closed. 


We remark that a rational number r is an integer if and only ifr € Zp 
for each prime p, and so this theorem may be regarded as a “local version” 
of Theorem 1.1. The proof shows that the local version follows at once from 
Theorem 1.1, and a simple bit of topology. 


Proof of Theorem 1.2. Let r = m/n, k, and (;) = s/t be as in the statement. 
Suppose a prime p divides ¢. If p does not divide n, thenr € Zp, so s/t € Zp, 
which is false. Thus each prime that divides ¢ divides n. 


4. POSTSCRIPT: AN “ELEMENTARY” PROOF 


A special case of Theorem 1.2 is that (7) is integral for all negative integers n and 
0<k EN. It is possible to prove this directly by extending Pascal's triangle 


backwards and using induction (we show a tilted version, in which the row 
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corresponding ton € Z has the entries (7), (k € N) beginning with (5) = 1, 
(7) =): 


1 —3 6 -10 15 
1-2 3 -4 5 
1 -1 1 -1 1 
1 0 0 0 0 
1 2 1 0 90 
1 3 3 1 0O 
1 4 6 4 1 


This case amounts to the following statement: 


Lemma 1. For each k € N, the product of any k consecutive positive integers is 
divisible by k!. 


Proof. When n € N, the numerator of (-"), in reduced form, is, up to sign, the 
product of the k consecutive integers starting at n. 


Given two integers a, b € N, we may factor b as b = rs, where r is relatively- 
prime to a, and s is a product of powers of primes that divide a. This fac- 
torization is unique: r may be specified as the largest factor of b such that 
gcd(a,r) = 1. We call r the part of b prime to a. 

A direct attack on Theorem 1.2 reduces it to showing the following gener- 
alization of Lemma 1: 


Theorem 4.1. Let d,k € Nand let r be the part of k! prime to d. Then r divides the 
product of the terms of each k-term arithmetic progression of integers having step d. 


Proof. For m € Z, let 
f(m) = m(m + d)---(m + (k — 1)d). 


We have to show that r| f(m) for each m € Z. 
Now gcd(r, d) = 1, so we may choose e € Z with de = 1 (mod r). Then 


e* f(m) =em(em + 1)---(em+k-—1) (modr). 


By Lemma 1, k! divides the right-hand-side, so r divides e* f(m), sor|f(m), 
since gcd(r,e) = 1. 
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THE MAGIC OF COMPLEX ANALYSIS 
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1. THE HISTORY OF COMPLEX NUMBERS 


Complex Analysis is all fun and games until someone loses an i. 


What are complex numbers? 


A complex number is a number which can be written in the form z = a + bi 
where a and b are real numbers and i is the imaginary unit satisfying i 2— 4], 
We call a the real part of z and b the imaginary part, written respectively as 


Re(z) = a and Im(z) = b. 


The set of all complex numbers is denoted C. The complex number a + bi can 
be visually represented as a vector (a, b) on a two-dimensional plane diagram 
called an Argand diagram. This diagram represents the complex plane and 
has two axes: the real axis and the imaginary axis. 


Im 


Figure 3.1: The complex plane represented on an Argand diagram. 


Where did it all come from? 


The very first mention of the notion of complex numbers goes back almost 
2000 years. In the first century, Heron of Alexandria was trying to calculate 
the volume of a pyramidal frustum. To solve it he had to find the square root 
of a negative number, and operation he deemed impossible, so he gave up. 
In the 1500s, some speculation about the square roots of negative numbers 
returned when formulas for solving 3°¢ and 4" degree polynomial equations 
were discovered. Mathematicians found that when using these formulae to 
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find real solutions of certain polynomial equations some intermediate work 
with negative square roots would be required. At the end of the calculation 
these negative square roots would cancel leaving a correct real solution. This 
problem led at least one mathematician to publish solutions to certain equa- 
tions without any clue as to how he actually found them! An example of this 
phenomenon, studied by the Italian mathematician Rafael Bombelli, can be 
seen with the equation 


x7 -15x-4=0. 


Employing the cubic formula gives a solution 
r= (2+ v-1) + (2- v-1) 


which formally involves imaginary terms. It simplifies to the valid real solu- 
tion x = 4. 

Another Italian mathematician Gerolamo Cardano worked with complex 
numbers around 1545. He attempted to find two numbers a and b satisfying 


a+b=10 and ab = 40. 
The solution to this problem is given by 
a=5+~¥7-15 and b=5-—~<V-15 


which we can easily find using the quadratic formula. These roots are certainly 
not real. It was becoming increasingly clear that imaginary numbers could not 
be avoided. 

The term “imaginary” itself was coined by René Desartes (not an Italian 
mathematician, but French!) in the 1630s to reflect his observation that “for 
every equation of degree n, we can imagine n roots which do not correspond 
to any real quantity”. This is an informal description of a key theorem in math- 
ematics, the Fundamental Theorem of Algebra, which more precisely states 
that every n" order polynomial with real coefficients has exactly n roots in C. 

Complex numbers gradually gained acceptance. John Wallis developed 
their geometric representation in 1685. The most rigorous development (and 
in asense justification) was given by the Irish mathematician Sir William Rowan 
Hamilton. He formally defined complex numbers as pairs of numbers (a, b) 
with associated operations 


(a,b) + (c,d) = (a+b,c +d) and (a,b)- (c,d) = (ac —bd,ac + bd). 


This purely algebraic construction sidesteps the philosophical issue of the mys- 
terious “number” i, which in Hamilton’s system is associated to the pair (0, 1). 
Many years later Hamilton extended this method of construction to a four di- 
mensional analogue of complex numbers, quaternions. 
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The Birth of Modern Day Complex Analysis 


Complex analysis is the branch of mathematical analysis that examines func- 
tions of complex numbers. Many complex functions are simple analogues of 
real functions, for instance 


f@)=2 and ge) =~, 


with the only difference being that we allow z to be complex. A number of 
prominent and well known mathematicians have contributed to complex anal- 
ysis, including Euler, Gauss, Riemann and, of course, Cauchy. 

There are several fundamental ideas in complex analysis. The derivative 
of a complex function is defined in a similar way to the real case, 


PO = tim LEA DELO 
h->0 h 

However it turns out that it is much harder for a complex function to be dif- 
ferentiable than a real function. The reason is that the limit above must be 
the same no matter which way h approaches 0. If h is real then it can only 
approach 0 from two directions: from the negative direction or the positive 
direction. However when h is complex it can approach 0 (the origin of the 
complex plane in figure 3.1) from uncountably many directions! 

This idea of it being difficult for a complex function to be differentiable is 
captured by one of the core concepts in complex analysis: the Cauchy-Riemann 
equations. Take some complex function f(z). The number z is complex and 
so can be written as z = x + yi for x and y real. The value of the function is 
also complex, and depends on x and y. We can thus write 


F(z) = FX + yi) = u(x, y) + 0%, yi 


for two functions u and v that each take in two real numbers and output a real 
number. 

Now suppose that f(z) is differentiable at a point z9 = xo + yoi. Then the 
functions u and v must be first order differentiable and must satisfy 


du 8 a a 
Oe ee (3.4) 


If this happens we have 
7) 

f'@o) = 5 + Si. 
ax 


The equations in (3.4) are the Cauchy-Riemann equations. If a function satis- 
fies them at zo then we say the function is analytic at zo. 

We said above that it is hard for a complex function to be differentiable, and 
the Cauchy-Riemann equations really demonstrate this. Suppose we create a 
complex function by randomly choosing two functions u and v — let’s make 
them infinitely differentiable to be easy on ourselves! — and define f(x + yi) = 
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u(x, y) + v(x, y)i. What are the chances that this randomly chosen function 
will satisfy (3.4) and hence be complex differentiable? It is almost impossible! 
For a given u there is only one associated v that will satisfy (3.4) everywhere, 
while there are more than uncountably many possible choices of v to pick from 
initially. So even though our two component functions u and v are infinitely 
differentiable our complex function is almost surely not even once differen- 
tiable in the complex sense! 

Now the question is, do we get anything in return for having such a strin- 
gent definition of differentiability? The answer is an emphatic yes. Let’s look 
at an important example, an important cornerstone of complex analysis called 
CIF — no, not the common household cleaning product, but the Cauchy Inte- 
gral formula. 

Let y be a simple closed curve in the complex plane, oriented counter- 
clockwise. Simple means the curve doesn’t intersect itself, closed that it has 
no endpoints, and the orientation just gives the direction we integrate over 
(recall for ordinary integrals we integrate left to right). So curve y might look 
something like this: 


Let f(z) be any complex function which is analytic at all points on and 
inside y. If Zo is inside the region bordered by y, then 


1. f © 


2ni Jy Z— Zo 


Ff (Zo) = dz. (3.5) 


So by just knowing the value of f(z) on y and the fact that it’s analytic we can 
find its value anywhere inside y. Even better, this formula extends to deriva- 
tives: for any n we have 


Ci ee F(Z) 
ete 201 [ (z aye 


The Cauchy Integral Formulae have many important consequences; we list 
three. 


o Ifafunction f is analytic at a given point, then its derivatives of all orders 
are analytic there too. 


o Let f be continuous on a domain D. If 


[ f(z)dz =0 


for every closed contour C in D, then f is analytic throughout D. 
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o Suppose that a function f is analytic inside and ona circle Cr, centered 
at Zo and of radius R. If Mr denotes the maximum value of | f(z)| on Cr, 
then 

n!M R 


If(zo)| = 


for every n € N. 


2. APPLICATION: REAL INTEGRALS 


Integrals are an integral part of mathematics, mind the pun. Actually comput- 
ing integrals is often important; however, most integrals cannot be solved by 
using “traditional” methods. Many of these real integrals (usually improper, 
meaning one or both of the limits of integration is +00) can be solved by using 
complex analytical methods. The primary tool are residues. 

To describe residues, let’s pick a complex function f that is analytic every- 
where except at certain points where it “blows up”. A good example is 


1 
f(z) = 7 


which is analytic everywhere except at z = 0. We say that f has a singularity 
atz = 0. 

A residue is just a complex number related to the singularity. First we pick 
any closed curve yo that encloses the singularity zo; for our example zo = 0 
and the unit circle will do as yo. Then the residue of f at the singularity zo is 
defined by 


Res(f, Zo) = a f(z)dz. 
Yo 


Residue Theorem 


Let y be a simple closed curve in C, oriented couterclockwise. If a function f 
is analytic on and inside y, expect possibly for a finite number of singularities 
at Z),...,Z,, then 


/ f(z)dz = 20i ° Res(f, zx). 
Y k=1 


Often this can give an easy way of evaluating real integrals which would be 
hard by conventional (real analytic) methods. Residues are unquestionably a 
useful tool for solving difficult real integrals, but sometimes the notation can 
be off-putting. Iam going to look at evaluating these integrals ina slightly dif- 
ferent way, using just the knowledge we already have and a little mathematical 
trickery. The best way to understand the method is with an example. 


Example 


We wish to evaluate 


cs dx 
=| G2 4+ D2 +4)" 
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This is the same as 


1% dx 
2 iz (x2 + 9)(x?2 + 4) 


Now we let ; 


~ PENG + A 
This function clearly has singularities at +27 and +37. Now, relabelling the 
integral, we see that we want to calculate 


f() 


: / - dz 
lim Se ae 

R>o J_R (z? + 9) (z2 + 4) 

The next step is to enclose the singularities above the real axis using the upper 
half circle as in figure 3.2. Let Cr be the semicircle of centre 0 and radius R > 3 
(chosen to larger than the modulus of the biggest singularity of f(z)). 


Im 
CR 
Re 
—R +R 
Figure 3.2: Upper Half Circle 
Now consider yr = [—R, R] UCR, oriented counter-clockwise. By the basic 


linearity of integration we have 


dz si dz dz 
yr (27 + 9)(z? + 4) = 2 (2? Oizo A) ai [. (z? + 9)(z? + 4) ee) 


We can calculate the integral on the far left using CIF, and we can also show 
that as R — oo the integral on the far right goes to 0. This will leave us with a 
value for the integral in the middle in the limit R — oo which is exactly what 
we want to calculate. 

We employ the Cauchy Integral formula on the left integral twice, once for 
each of the singularities enclosed by yr. How do we do this? Observe that we 
can factorize f(z) as 


1 
Zz — 2i)(z + 2i)(z — 3i)(z + 31) 


A aay 
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We are interested in calculating the integral of f along a small curve a that 
only contains the singularity z = 21. Now inside this curve the function 


1 
(z + 2i)(z — 37)(z + 37) 


g(z) = f(Z)@ — 2i) = 


is analytic as all of the remaining singularites are outside the curve. Now by 
CIF (3.5) we have, with zo = 2i, 


eai=— f 2 


- -dz, 
2mi Jy Z—2i 


But g(z)/(z — 27) is just f(z), giving 
g(2i) x 2ni = 1 f(z)dz, 


and so 


1 — 
ees (a; + 21) Qi — 312i + sm) SAG 


Similarly, the integral of a small curve around the singularity at z = 37 is—/15 
(work it out!) and this finally gives 


dz ee: a 
ye (22 +9)(22 +4) 10 15 30° 


Now we want to show that the limit as R — oo of the right integral is zero; 


; dz 
lim / —_________ = 0 
R>00 JCp (Z7 + 9)(z? + 4) 


This is a basic exercise in bounds: 


/ dz </ |dz| 

eg PUG 4) Jeg e  llZe | 
<f \dz| 
~ Jer [1271 — 91127 |All 


1 
~ (R2 = 9)(R2 — 4) [. foe 


wR 
= (R? —9(R2—4) —>0,as R> ow. 


Subbing all this into (3.6) we get 


rf a +0,as R > co 
30. Jer (22 4+9)(22 +4)” 


and hence 


r=[- dz La 
Jo (22 49)(22 +4) 230 60° 
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We have solved the original integration problem. 

This almost ad-hoc method for solving these types of integrals is very ef- 
fective and relatively straightforward. It requires no more knowledge than we 
have already acquired and the simplicity of breaking the integral we want to 
find into integrals we can find a solution to makes this method very attractive. 

The idea that a real integral can be solved by using complex analysis meth- 
ods is quite amazing. Complex analysis has a vast amount of applications not 
only in mathematics and applied mathematics but also in physics and engi- 
neering. 


3. OTHER APPLICATIONS 


There are many applications of complex numbers and complex analysis. In 
electrical engineering, complex numbers are perfect for treating amplitudes 
and phases correctly as AC currents and voltages are combined. They are 
also used in signal processing, which has applications to telecommunications, 
radar (which assists the navigation of aircraft), and even biology (in the analy- 
sis of firing events from neurons in the brain). They are a fundamental aspect 
of the physical field of quantum mechanics. Imaginary (complex) numbers 
help form the descriptions of electronic states in materials which lead to ap- 
plications in optics. Complex numbers and complex functions have an impor- 
tant role in the field of fluid mechanics. Potential flows can be dealt with very 
easily through functions of complex numbers. 

In number theory, Riemann introduced revolutionary ideas into the sub- 
ject, the chief of them being that the distribution of prime numbers is inti- 
mately connected with the zeros of the analytically extended Riemann zeta 
function of a complex variable. In particular, Riemann introduced the idea of 
applying methods of complex analysis to the study of the the prime number 
counting function (x), which gives the number of primes less than or equal 
to x. 


4. CONCLUSION 


There are many applications, some unexpected, of complex numbers and anal- 
ysis in the real world. It has come from rocky origins and skeptical develop- 
ment, to its initial mainstream acceptance and now worldwide appreciation. 
Many notable mathematicians have contributed to this intriguing branch of 
mathematics, but it doesn’t stop there. The future of complex analysis is very 
much active. There is a lot of current research in this area, especially in Hy- 
percomplex analysis and Quaternionic analysis. In fact, if you fancy winning 
yourself one million US dollars, see if you can solve one of the millennium 
problems, the Riemann hypothesis. 


Question 


The distribution of the prime numbers among all natural numbers does not 
follow any regular pattern; however, Riemann observed that the frequency of 
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prime numbers is very closely related to the behavior of an elaborate function 


1 1 1 
called the Riemann Zeta function. The Riemann hypothesis asserts that all 
interesting solutions of the equation 


S(s) =0 


occur when the real part of s is precisely one half. This has been checked 
for the first one and a half billion solutions. A proof that it is true for every 
interesting solution would shed light on many of the mysteries surrounding 
the distribution of prime numbers. 

Did you know i! is a real number? Check it out! 
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THE FALSE POSITIVE PARADOX 


MICHAEL HANRAHAN 
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The false positive paradox is a counter-intuitive statistical result whereby a 
positive test result is more likely to be a false positive than a true positive. 
This occurs when the incidence of the thing you are testing for is lower than 
the false positive rate of the test. This phenomenon can cause confusion when 
dealing with medical test results for very rare conditions. The result of it is 
that when you get a positive result for a rare condition the actually probability 
that you do have the condition is still quite small even for very accurate tests. 

There are screening tests for many medical conditions and diseases that 
doctors use to diagnose people who may be at risk. For very rare diseases the 
vast majority of people screened will get a negative result and generally the 
risk of getting a false negative result is low. However when someone gets a 
positive result it does not necessarily mean that the patient actually does have 
the disease or condition. This causes a lot of anxiety, stress and additional 
financial costs for further testing. On top of this, unnecessary additional test- 
ing can be painful and lead to complications. For instance, a prostate biopsy 
conducted after a positive PSA test can have complications such as infection, 
bleeding into the urethra or bladder, inability to urinate, bleeding into the rec- 
tum, or allergic reactions to the anaesthetic. 


1. BAYES’ THEOREM 


In order to calculate the probability of actually having a true positive result 
when the test is positive statisticians use Bayes’ theorem. To use Bayes’ theo- 
rem in relation to medical screening the following information is needed: 


o The prevalence of the condition/ disease in the population being tested. 


o The sensitivity (actual positives correctly identified) and specificity (ac- 
tual negatives correctly identified) of the screening test. 


Mathematically, the theorem is expressed as: 


P(B|A) P(A) 


where the event A signifies having the condition, the event B signifies having 
a positive test result, and hence 


o P(A|B) represents the chance of having the condition given a positive 
test result. 


o P(B|A) represents the chance of a true positive; i-e., the positive test re- 
sult being correct. 
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Figure 4.1: A represents the whole population to be tested, B represents the 
number of people who have the condition and C represents the people who 
get a positive test result. So it’s clear that most of the people who do have the 
condition do get a positive result, but so do a lot of healthy people. The false 
positive paradox arises because C is at least twice as big as B(\C. 


o P(A) represents the incidence of the condition in the population, 
o P(—A) represents the chance of not having the condition (i.e. 1 — P(A)) 


o P(B|—A) represents the chance of a false positive; i-e., getting a positive 
test result but not actually have the condition. 


In words the theorem says 


number of true positives 


chance that a positive result is correct = es 
number of total positives 


When using the theorem it can be helpful to write the variables in table format. 


Has condition, P(A) | Doesn't have condition, P(—A) 
Positive Test Result | True Positive, P(B|A) False Positive, P(B| — A) 
P(B) (test sensitivity) 
Negative Test Result False Negative True Negative 
P(-B) (test specificity) 


The false positive paradox arises because medical screening tests are never 
100% specific (a screening that was specific would never give false positive) 
and for rare diseases a small percentage of a large number of people (number 
of false positives) can be larger than a high percentage of a small number of 
people (number of true positives). See the visualization in figure 4.1. 

We now look at some real-life applications of Bayes’ theorem to common 
screening tests. 
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2. FIRST EXAMPLE: BREAST CANCER 


The prevalence of breast cancer in Ireland is approximately 0.2% in women 
aged over 40. The screening mammography has a sensitivity of 69% and a 
specificity of 94%. This gives rise to following the table. 


Has Breast Cancer (0.2%) | No Breast Cancer (99.8%) 


Positive Test Result 69% 6% 
Negative Test Result 31% 94% 
We then get 


True Positive = 69%«0.2% = 0.00138, 
False Positive = 6%*99.8% = 0.05988. 


We then have that the probability of having Breast Cancer after receiving a 
positive test result is 


true positives 


= 2.25%. 
number of all positives 


This is a shocking result and is the source of much controversy as to whether 
mammograms are even worth doing. Approximately 977 people out of 1000 
people who have had a positive mammogram result do not actually have can- 
cer. About 10% of those people will be sent for additional testing. 14% of that 
10% will be sent for a biopsy, and of this group 65% will still not have breast 
cancer owing to the fact that the biopsy test has its own false positive rate. 

Another result of counter-intuitive probability states that over a decade of 
annual mammogram screening about 46% of women will receive a positive 
result, the vast majority of which will be false positives. 

We have now clearly shown some disadvantages of mammogram screen- 
ing in terms of unnecessary costs and anxiety for patients. The question that 
remains is why do we tolerate such high false positive rates? Surely this test 
can be refined to give fewer false positives? 


3. ROC CURVES 


The reason why the test performs to these standards can be explained by the 
Receiver Operator Characteristic (ROC) curve of the test. The curve demon- 
strates that there must be a “trade-off” between sensitivity and specificity. This 
arises because, in general, there is no perfect marker for a disease. The levels 
of the marker found in any individual with or without the disease will follow 
a general pattern of distribution. The problem is that overlapping of the two 
distribution curves of healthy people and people with the disease occurs so 
that at a certain concentration of the marker it is difficult to tell whether or not 
the person has the disease. 
And example distribution is illustrated in figure 4.2. The regions are: 
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Threshhold 


Healthy people 


People with disease 


Number of people 


Level of test marker 


Figure 4.2: Distribution curves for a particular disease. 


1,3. Healthy people who receive true negative test results. 
2,4. Healthy people who receive false positive test results. 

3. People with the disease that receive false negative test results. 
4,5. People with the disease that receive true positive test results. 


The threshold is an arbitrary cut-off point designated by the creator of the 
test above which any level of marker detected will yield a positive test result. 
The ideal threshold point maximises the number of true positives and min- 
imises the number of false positives. However in some instances it may be 
more useful to have a high specificity. 

In order to visualize the effect of moving the threshold point (i.e., the trade- 
off between sensitivity and specificity) we use the ROC curve (figure 4.3). 


4. SECOND EXAMPLE: HIV 


Human immunodeficiency virus (HIV) is a virus that infects the human im- 
mune system and leaves the infected person susceptible to life-threatening op- 
portunistic infections and cancers. It can cause acquired immunodeficiency 
syndrome (AIDS). 

The prevalence rate of HIV in adults (aged 15 to 49) in Ireland is 0.2% (2009 
estimate from UNICEF). An oral home-use HIV test can be used to detect 
whether or not HIV antibodies are present in someone who is at risk of be- 
ing infected with the HIV virus. For this test, 99.9% of negative results are 
true negatives (i-e., the person does not have HIV antibodies) while 91.7% of 
positive results are true positives. We put these numbers in table form: 
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True positive rate (sensitivity) 


False positive rate (100-specificity) 


Figure 4.3: ROC curve: the more sensitive we want the test to be, the less 
specific it will be. 


HIV Positive (0.2%) | HIV Negative (99.8%) 


Positive test result 91.70% 0.10% 


Negative test result 8.30% 99.90% 


What is the likelihood that you have HIV if you get a positive result? As 
before we calculate the numbers: 


true positive = 91.7%*0.2% = 0.001834, 
false positive = 0.1%*99.8% = 0.000998. 


We then get that the probability of being HIV positive is 


, is 
SEE POSEYES — 0.001834 /0.002832 = 64.75% 


all positives 


Therefore, there is only a 64.75% chance that you truly are infected with the 
HIV virus if you get a positive result. In other words, roughly 35 people out of 
100 who get a positive result using this test will actually not be infected with 
HIV. This is a counter-intuitive result as the test originally appeared to have 
a 91.7% accuracy rate. If someone does get a positive test result it cannot be 
inferred that they do have HIV and further testing is required. However it is 
easy to see how this could cause a lot of unnecessary anxiety and stress on a 
person. 

Clearly this home test was designed to give a true negative result to peo- 
ple who do not have HIV. Using the same method as above it’s found that a 
negative test result means there is a 99.99% chance that you don’t have HIV. 
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So this wasn’t a true example of the false positive paradox but it’s still a 
very interesting example and highlights the importance of not over-reacting 
to a positive test result. 


CONSTRUCTING THE REAL NUMBERS 


DARON ANDERSON 
Trinity College Dublin 


Our goal is simple: to show that the set of real numbers, which we use every 
day and which seem self-evident from the properties of observable space, can 
be defined in a purely abstract way that is free from self-contradiction. 

The real numbers can be intuitively understood in terms of continuous 
amounts, or as positions on the number line. But if we try to build a definition 
around this concept, removing all references to an external world, it is hard to 
see what are we left with. What is an amount if not an amount of something? 
The amounts on their own, considered without the associated operations of 
addition, multiplication, order, and absolute value, are of little interest. In- 
deed, in this form, the only notable property is the set’s cardinality, a measure 
of how large it is. 

Thinking of the set R like this, with no way of comparing or manipulating 
the elements, gets us nowhere. The characteristics of the real numbers that 
we are interested in - those we would like to prove make sense - include the 
manner in which they interact under the rules of arithmetic, how these rules 
interact with each other, and the specific way in which we have the set inclu- 
sions N Cc ZC Q CR. These well known simpler subsets of the reals can be 
defined in a purely symbolic manner. We will do this, if only to illustrate the 
comparative ease with which it can be done. 

The natural numbers can be defined rigorously as elements of the form 
1+1+4...+ 1, with addition defined in the obvious manner. We define mul- 
tiplication by setting 1-1 = 1, and by imposing the distributive law: that 
a-(b+c) =a-b+a-c, for every a,b,c €N. This extends to Z simply enough 
by introducing one new element —1 such that 1 + (—1) = 0 and again asserting 
that multiplication distributes over addition. 

We then construct the rational numbers Q as elements of the form § where 
a and b are integers and 5 is nonzero. We then define 


5 = G tomean ad = be, 
a c ad +cb 
eee me (tie 2 
pt q to be equa to ‘dc 
2 te bee ual to “= 
bd d bd’ 


and we're done. 

Attempts to extend this technique to R provoke some immediate questions. 
For example, what exactly does it mean to say that 2.71828... has an infinite 
number of decimal places? This question is especially difficult since 2.71828... 
is not so much a mathematical object as a method of representing one. 

We could conceivably define the reals in an algorithmic way, as maps of 
the form x : N — {0,1,2,3,4,5,6, 7,8, 9} which assign to each n what we 
understand to be the n decimal digit of x, along with rules for adding and 
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multiplying such maps. However, unlike in the case of rational numbers, these 
algorithms must be infinite. And it is difficult to define the result of an algo- 
rithm which never terminates. We can sidestep this problem by giving, instead 
of a rule for calculating the entire decimal expansion of x + y, say, a recursive 
procedure which will produce any chosen (x + y)(7) after a finite number of 
iterations. This is a valid yet unappealing definition of the real numbers.* 

Overall, there are three problems with this approach. The first is that it 
places undue importance on the decimal representation of numbers: math- 
ematically speaking, there is nothing special about this way of writing real 
numbers. We could just as easily write the number in binary, or indeed any 
integer or non integer base, or use continued fractions. The second problem is 
that this system does not capture the fact that the decimal expansion of a real 
number is not necessarily unique. For instance 


3 1 
Las 5-15 3-5 = 363.333...) = 0.999... 


and therefore 1 = 0.999... Finally, we want to be able to prove things about 
the real numbers and use them to study other concepts without always having 
to refer back to the nitty-gritty of how addition and multiplication work. We 
want to be able to say let x be a real number without worrying about how to 
write x down. In order to do this, we would most likely have to prove that 
this construction of R is equivalent to some other more natural one. This then 
renders the original construction redundant. 


1. FIRST CONSTRUCTION: AS A COMPLETION 


What more natural constructions exist? The first one here relies on the fact 
that for each real number, there are rational numbers arbitrarily close to it. We 
must make rigorous the concept of estimating reals by rationals. 

The first step is to present the rational numbers differently. For each g € Q 
let g denote the rational sequence (q,q,q,...). We can define addition and 
multiplication pointwise on the space of these sequences. This gives us a struc- 
ture which behaves identically to the rational numbers. 

What happens if we allow non constant sequences? Perhaps if we allow 
all sequences we get a collection of objects that behaves like the real numbers? 
We can certainly allow non constant sequences, but it turns out the space of all 
sequences is simply too big to be R. For instance there is only one element of 
R that cannot be inverted - namely 0 - while there are many sequences other 
than 0 which are non-invertible. For instance, there is nothing which when 
multiplied by (0, 1, 1,...) yields 1. But we want R to contain only one element 
0 with this non-invertibility property. 


*TEditor’s note]: This method of construction is also in a certain sense impossible. Just as the 
set of all numbers with finite decimal expansion is a countable set, so is the set of all “finite algo- 
rithms”; that is, algorithms that can be defined using a finite number of instructions. As the set 
of real numbers is uncountable it follows that almost all real numbers cannot be produced using 
a finite algorithm. These numbers are termed uncomputable. For a popular account of this issue 
see the second chapter of Roger Penrose’s book The Emperor's New Mind. 


36 CONSTRUCTING THE REAL NUMBERS [Summer 


However, there is a special category of sequences which suits our purpose. 
We allow exactly those sequences which tend towards what we would naively 
regard as a real number. These are called Cauchy sequences. For example 
(3, 3.1, 3.14, 3.141, 3.1415, 3.14159, ...) is allowed. The terms in this sequence 
get arbitrarily close to each other, and we would intuitively expect this means 
the sequence converges. It is true that over the real numbers, this sequence 
converges. But since we have not defined the real numbers yet, such a state- 
ment is so far meaningless. 

Clearly more than one sequence can approximate the same element of R 
and so we must consider some of the allowed sequences to be the same. We 


would like to allow (xm)?°_, and (m)?°_, to be equivalent if 


lim Xn = lim ym, 
m—->oo m—->Co 


but since these limits may not be defined the best we can hope for is that 
lim (%m — Ym) = 0. 
m—>Co 


This makes sense since every Xm — Ym is a rational number. 

If a number x is rational, we can represent it in this space by x. If we want 
the number to be irrational, we can use the sequence (%m)%_, which has x, 
as the decimal expansion of x truncated to n decimal places. This is clearly a 
rational sequence. 

In order to suggest the rules of arithmetic on these Cauchy sequences work 
as intended, we will assume that R is already defined and that, rather than 
acting as a definition, this construction simply mimics it. Now x is represented 
by (Xm)p,—, With limit x, and y is represented by (ym)?°_, and with limit y. 
Then x + y is represented by sequence (Xm + Vm)°°_, as 


lim (Xm + ym) = lim Xm + lim ym = x+y. 
m—>oo m—>oo m-—>oo 


The rule for multiplication is defined similarly. We finally need to define the 
ordering on the reals. We define x < y to mean that there exists M € N such 
that 

Xm —Ym < Oforallm > M, 


which allows us to order these sequences. 

What we have done here is encoded the information required to represent a 
real number in a sequence, and then used the sequences themselves as our ba- 
sic objects. Thus we have translated the problem about an infinite decimal ex- 
pansion into one about infinite sequences, which are already well understood. 
This approach can be generalized to form what is called the completion of a 
metric space. You might hear the real numbers described as the completion of 
the rationals. This abbreviates the entire construction above. 


2. SECOND CONSTRUCTION: DEDEKIND CUTS 


While our previous construction was topological, the second involves more set 
theory. Dedekind cuts are named after their inventor, illustrating the signifi- 
cance of the individual contribution to this general approach. 
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Dedekind noticed that a real number can be completely specified by giving 
two sets: the set A of numbers which are less than or equal to it, and the set B 
of those which are strictly larger than it. Furthermore, he noticed that nothing 
is lost if A and B only contain the relevant rationals. He then flipped this 
realization on its head and defined a real number as a pair of sets of rational 
numbers (A, B) with the following properties: 


1. Aand B are disjoint. 

2. Between them, A and B exclude at most one rational number. 

3. B does not have a minimum element. 

4. If A or B contains q; and q2 and gq; <q < q2 then it contains q as well. 


What should always be remembered is that when defining a Dedekind cut 
we never assume there is some point lying between A and B on the rational 
number line. We use the cut as an object in its own right. 

The representation of an element g € Q as the Dedekind cut (A, B) is par- 
ticularly simple: we set 


A={xeQ:x<gsandB={xeQ:x>q}. 


Thus if A has a maximum, the cut (A, B) represents a rational number. Let’s 
try something harder. How about /3? We can set 


A= {x €Q:x <0V x? <3} and B= {x €Q:x SOA x? > 3}. 


where v denotes logical or and A denotes logical aNp. Here there is no confu- 
sion over which representations are equivalent. We simply say that (A, B) = 
(C, D) if and only if A = C and B = D as sets. Then we define the operations 
of addition and multiplication as follows: (A1, B1) + (A2, Bz) = (A142. Bi+2) 
for 


Aj42 = {a1 + a2: a, € A, A a2 € Az} and 
By42 = {b, + b2: by € By A bz € Bo}. 


For (A, B) rational, we define —(A, B) in the obvious way. Otherwise we set 
—(A, B) = (—B,—A) for — A = {-—a: a € A}, etc. 


This permits subtraction. 
For ordering we denote by (A1, B1) < (Az, Bz) that Ay © Az. This allows 
us to define the sign of a cut, 


1 if0 <(A,B)A(A,B) £0, 


sen(A,B)= 4 —-1 if(A,B) <0A(A,B) £0, 
0 if(A,B)=0, 
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which in turn allows us to talk about positive and negative numbers. 
Multiplication is first defined for positive numbers through 


(Ay, By): (A2, Bo) = (A1-2, B12) 
where 


Ay2 = {41-2 : a, € Ay A a2 € Ap Aaj, a2 > O} U {x € Q: x < O} and 
By2 = {bi - bo: by € By Abo € Bo}. 


Multiplication is extended to all numbers by setting 


(Ai, Bi) - (Ao, B2) = sgn(A1, Bi) - sgn(A2z, Bo): 
(sgn(A1, Bi)(A1, Bi) - (sgn(A2, B2)(A2, Bo). 


These definitions give all the expected properties of the real numbers, such 
as commutativity and distribution of multiplication over addition. Thus we 
have built up an algebraic system which interacts with the rational numbers 
in the required way. 

Now it is not immediately obvious that under the “continuous amounts” 
interpretation, that each real number can be represented by a Dedekind cut. 
Certainly those with nice properties such as x? = 2 can be, but do all real 
numbers have such clean rules for deciding exactly where they live on the 
continuum? Let’s go back for a moment to our definition of R in terms of 
Cauchy sequences. It should be obvious that for each real x, some rationals are 
less than or equal to x, while some are greater. Moreover these two sets satisfy 
the Dedekind cut axioms. But describing them without explicit reference to 
x, that is the challenge. Another possible criticism of Dedekind’s approach is 
that it doesn’t say what a real number looks like, or how to write one down 
on a page. This is a blessing in disguise, since it leaves us free to use whatever 
notation we like, provided it makes logical sense in this framework. 


3. THIRD CONSTRUCTION: ALGEBRAICALLY 


The next construction of R requires some background, and is if anything more 
abstract than the previous. But since the following required definitions are 
already ubiquitous across mathematics, if you have not seen them before you 
will again soon. 

The first is the notion of a group, which is a generalization of the rational 
numbers. It is a set of elements which can be “multiplied” together in a man- 
ner which satisfies many familiar properties. Common examples of groups 
include Z and Q under addition, and Z \ {0} and Q \ {0} under multiplication. 
(We must omit 0 in the last examples as the group operation is required to 
be invertible.) Slightly less common examples are the set of invertible square 
matrices of a given size, or injective maps from a set onto itself under compo- 
sition. 

Formally, a group is defined as a pair (G,-) consisting of a set G and a 
binary operation - on G fulfilling the following axioms 
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1. Foralla,b € G,a-b € G. (Closed under multiplication.) 
2. For all a,b,c € G, (a-b)-c =a-(b-c) (Multiplication is associative.) 


3. There exists e € G such that for alla ¢ G,e-a =a = a-e (Identity 
element.) 

4. For alla € G there exists a~! € G such thata7!-a =e =a-at. 

(Invertibility.) 


If, in addition, a-b = b-a for every pair of elements a,b € G then G is called 
an Abelian group and we generally write + instead of -. 

A field (F, +, -) consists of a set F equipped with two binary operations + 
and - , where (F,+) is an abelian group with identity element denoted by 0, 
and (F \ {0}, -) is an abelian group with identity element denoted by 1, and 
the field elements also satisfy the distributive law a-(b +c) =a-b+a-c. 

A field is said to be totally ordered if it is equipped with an abstract binary 
relation, denoted <, such that 


1. Foralla,be F,a<bvb <a. (Totality.) 

2. Foralla,be F,a<bAb<a>a=b. (Antisymmetry.) 

3. Foralla,b,c¢ F,a<bAb<c >a <c. (Transitivity.) 

4. Foralla,b,c€ F,a<b3>a+c<b+c. (Order commutes with sum.) 


5. Foralla,b € F,0 <aA0 <b> 0 < a-b (Order commutes with 
product.) 


Since we can add the | of the field, or its negative, to itself n times, we can gen- 
erate a copy of the integers living inside our ordered field. Phrased precisely, 
each ordered field contains a subfield homomorphic to Z. 

Finally, we define an Archimidean ordered field: an ordered field is said to 
be Archimidean if, for any element x of F, there is a smallest integer greater 
than x. 

The point of all of these definitions is that it can be shown there exists a 
maximal Archimidean ordered field, meaning one which contains a copy of 
any other Archimidean ordered field. This unique field can be used as a defi- 
nition of the real numbers. 

This approach is based solely on abstract algebra. Most textbooks on real 
analysis will begin with a list of axioms for the real numbers. And if you look 
hard enough at the definitions above, you can recover most of these axioms. 
What the algebra gives us is that these axioms define an unique structure. This 
is not just our favourite Archimidean ordered field; it is the only Archimidean 
ordered field! 
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4. FOURTH CONSTRUCTION: SURREAL NUMBERS 


Our final construction is called the Theory of Surreal Numbers. It differs from 
the previous three in two ways. First, it does not contain any undefined notion 
of equality. Second, this schema does not rely on us already having established 
the natural numbers. It simultaneously generates the naturals, integers, ratio- 
nals, and reals as well as giving a rigorous definition of infinite numbers and 
infinitessimals. 

This system deals with the field of numbers which are defined recursively 
as follows: (@|@) is a number and so is any symbol of the form (X|Y) where 
X and Y are sets of numbers such that if x € X and y € Y thenx ¢ y. 
This becomes more complicated when you realise we haven't yet defined what 
the statement x < y means. This total ordering of the set of numbers is also 
defined recursively. We first let (0|@) < (@|@). Then for x = (Xz|Xpr) and 
y = (Y¥z|Yr) we write x < y if every x, € Xz satisfies y < xz and every 
yr € Yp satisfies yr £< x. But enough theory. Let’s actually construct some 
numbers. First, the naturals. 


o Let 0 abbreviate (6|9). 

o Let 1 abbreviate (0/0) = ((@|®)|@). 

o Let 2 abbreviate (1|@) = (((@|)|@)|@). 

Oi. 

o Let the numeral 1 + 1 abbreviate the number (n|9). 
Next the negative integers. 

o Let -1 abbreviate (@|0) = (@|(@|@)) 

o Let -2 abbreviate (@| — 1) = (8(G|(@|))) 

Or. 

o Let the numeral —(7 + 1) abbreviate the number (@| — 7) 


Note that we really should write ({0}|@) instead of (0|%) and so on, but we will 
omit the curly braces for the sake of convenience. 

Before we go any further, let’s check these satisfy the order relations we are 
looking for. Returning to our definitions, we shall prove 0 < 1; that is, 


0 = (@|9) = (Xz, |Xr) < (Yz|Yr) = (0|) = 1. 


Indeed, since Xz, and Yr are both empty, both properties are vacuously true. 
As that was rather easy, we will try | < 2; that is, 


1 = (0|9) = (Xz|Xr) < (Y/Y) = (1|9) = 2. 


Now Yz is still empty, so the second property is true. For the first we need to 
prove that for every xz € {0}, (1|@)£xz. Since the only element of {0} is 0, we 
then need to prove 0 < 2. The proof of this is vacuous as for 0 < 1. 
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Using standard induction we can establish the integers and show their or- 
der behaves as desired. Let’s go further and define powers of 4. 


o Let 5 abbreviate (0|1) 

o Let 4 abbreviate (0|4) 
©... 

o Let eas abbreviate (0| 47) 


Now we need a way to add and multiply things like (X|Y). Again this is 
done recursively. First we set 


0+ x = (O|9) + (XL|Xr) = (XL|Xr) = x. 
We then define 
X+Y = (Xz|Xp) + (Yi|Yr) = (Xi + y,x + Yt|Xr + y,x + Yr) 
and allow subtraction through 
—0=0 and — X = —(Xz|Xr) = (-Xz|— XR). 
We define multiplication by setting 
XY = (X1|Xp) + (Yz|Yr) = (Z1|Zr) 
where 


Zi ={y- Xt +tx-Y¥,—-Xi-Yi,y-XrR+x-Yr—Xr-Yr}and 
Zr={y-Xpt+x-Yr—Xzi-Yrex-Yrt+y-Xr—XpR-Yz}. 


These rules of arithmetic are a little difficult to follow so, if you like, it may 
help to try to add or multiply some small numbers. 

We can use the integers and powers of a half to construct what are called 
the dyadic fractions, fractions which have a power of two as their denominator. 
We can then define a real number in a way similar to a Dedekind cut, choosing 
r = (X|Y) where X is the infinite set of dyadic fractions we want to be less 
than or equal to r, and Y the set of dyadic fractions we want to be greater. This 
determines r. 

Once again, not every (X|Y) symbol represents a unique number. What we 
can say is that if (Xz|Xr) < (Yx|Yr) and (Yz|Yr) < (Xz|Xpr) then it makes 
sense to consider (X;|XR) and (Y,|Yr) as equivalent. Since we haven't used 
the notation yet, we can denote this relation by (X,|Xr) = (Yt|Yr). 

Remember that we originally defined (X|Y) by induction. It is clear that 
if we have an infinite list of numbers which we wish to construct, then we 
should be able to do this through some form of induction involving the (.|.) 
operator. But one property of the real numbers is that they are uncountable 
meaning they cannot be arranged on a list. Considering this, it seems unlikely 
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that we will be able to build all of R in this manner using standard induction. 
For this we will require an axiom governing transfinite induction. This allows 
chain-of-effect type proofs where our dominoes are not necessarily discrete. It 
allows such deductions as that (1, 2, ...|@) is a number, even though it can not 
be reduced to something in terms (0|@), as every integer can. 

While the soundness of transfinite induction seems like a strange thing to 
believe, it is an equivalent assumption to about a million other things, some of 
which seem more intuitive, including the notorious Axiom of Choice. 

The theory of Surreal Numbers does not rely on any previously constructed 
mathematical objects, though it does require some logical axioms to function. 
And in order to follow the proofs, one requires explicit knowledge of set the- 
oretic axioms, which explains why it is not often taught to undergraduates. 
This is not to say that the previous three constructions do not contain set the- 
ory; rather in those cases it is hidden below the surface, so that it seems more 
like common sense than anything else. 


5. AND ON... 


The real numbers can be defined in many equivalent, but seemingly unrelated 
ways. They may not be an inherent property of the natural world. But the 
fact that there are so many different coloured roads to the same destination 
suggests that they live somewhere at the heart of our mathematical logic. 

So once more from the top: let x be a real number... 


ALL KINDS OF PI 


JAMES FENNELL 
University College Cork 


As a reader of this article you will be familiar with the constant z. In fact, you 
are probably so familiar with z that you've perhaps become a little complacent 
about it. You know that z is irrational, but — unlike /2 — you're not exactly 
sure how you'd go about proving it. You know that, as an irrational number, z 
has an infinite and non-repeating decimal expansion, but will likely struggle to 
give anything more than the first 3 or 4 digits. But you’re at least sure that your 
personal inability to recall these specific facts is irrelevant: of all mathematical 
concepts z is surely the one with the most static and objective existence, its 
eternal character always there in the background ready to be recalled when 
needed. The goal of this article is to show that, in thinking that, you are in a 
specific sense wrong. 

In the following pages we will see that, rather than being some universal 
unchanging entity, the idea of 7 as we know it is inextricably linked with our 
personal notions about the geometry of space — and is hence subject to same 
kind of arbitrariness that those notions are. We will realize that statements 
like 


15 
z=3 ort =v15 tag (6.7) 


are not necessarily the ramblings of numerically illiterate engineers (though 
you should always double-check) but instead represent a different and more 
general way of looking at the 2 concept. Indeed, we will see that (6.7) are all 
legitimate values of a but that — and I say this now to avoid the impression 
that things will soon become a free-for-all — we reasonably cannot have 


21 
x= ~/8 or sea or 


Read on to see why! 

It is natural to ask if this tinkering with such an age-old concept is anything 
more than mathematical trouble-making. It will actually turn out that, apart 
from being intriguing for its own sake, the idea of generalizing 2 will leave us 
with a charming if limited mathematical tool. 


1. FROM WHERE COMES JT? 


To begin, we need to think about where z comes from. There are multiple 
ways of defining it. The first geometric way is to reason that, as a two-dimensional 
plane figure, the area of a circle should be proportional to the square of its ra- 
dius, and set z to be the proportionality constant to give 


area of circle of radius r = mr’. 
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We can use integration theory to give precise meaning to the term area and, 
indeed, us as the definition itself 


1 
mina f V1—x? dx. 
0 


This last formulation is nice because it does not depend on intuitive geometric 
notions which, as seen frequently in the history of mathematics, can often be a 
barrier to mathematical rigour. The integral definition can be made purely for- 
mal. Keeping this mindset we can design another definition free of geometry: 
formally define the function 


x? x4 it x 
cos(x) = 1- — + ——++-= > C1" 
n=0 


2n 


21° 4! (2n)!’ 


which makes sense as the infinite series converges for every x € R by the 
ratio test. Then set z to be twice the smallest positive number x satisfying 
cos(x) = 0. This definition is interesting because there is no a priori reason to 
believe that this equation actually has a solution and hence that the number z 
exists. 

However when people are thinking about z it is rarely in terms of the zeros 
of the cosine function. The most common way of viewing the origin of z is in 
its relation to the circumference of a circle. We can define 


T= [circumference of any circle of radius r]. (6.8) 
For this to make sense we need to define what we mean by “circle of radius 
r” and “circumference”. When we do this we become aware of the certain 
arbitrariness in our value of x. Let us see how. By a circle centered at x of 
radius r we mean the set of all points a distance r from x, where by the distance 
from points x = (x1,.x2) and y = (yi, y2) in R? we mean 


Ix — yllp = VG —y)? + Gr — yn). (6.9) 


But why should we be obliged to use this specific notion of distance? 

The best known alternative to the standard distance is the taxicab or Man- 
hattan norm, where we define the distance between the two points x and y in 
the plane as 

l|x — yl], = [xa — yal + [x2 — yal. (6.10) 


The reason for the subscript 1 will become apparent later on. Now this alter- 
native distance has all of the usual properties that we should expect from any 
kind of distance function. It is always positive and it is zero if and only if x 
and y are the exactly same point. It satisfies the triangle inequality: if z is a 
third point then getting from x to y by passing through z will always take at 
least longer than going from x to y directly. We can write this symbolically as 


lx yh Sle —2hh + lle - yh. 
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Though (6.9) has the strength of the Pythagorean theorem behind it, there 
are scenarios in which (6.10) is more appropriate; namely if you are walking 
through Manhattan and can’t travel “as the crow flies”. 


@ 


@ @ 
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& & é 
x & > & s & 
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Pe nS ie eee rae ee e = 42" street 
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Sa SSS aie aS ste ees - ryan Public 41° street 
Park -——— 
Library 


40" street 


we----- 


lx — yl], = |7—5| + |40 — 41| = 4 streets. 


If we decide to amend our notion of distance, what will the “circles” look 
like? To avoid confusion when working in the more general setting, we use the 
term ball instead of disc, and the term boundary of the ball instead of circle. 
The boundary of the unit ball —- the ball centered at the origin and of radius 1 
— is defined analogously to the usual case as the set 


{x € R?: ||x —O|| = ||x]| = ee 
If we denote the boundary of the unit ball with respect to || - ||, by Bi, then 
By = {x = (x1, X2) € R? : |x| + |x2| = 1}. 


In the upper right quadrant of the plane, x; > 0 and x2 > 0, so points on the 
unit ball will satisfy x2 = 1—x,. In the upper left quadrant we have x; < 0 and 
X2 = 0, and then as |x;| = —x, points will satisfy x. = 1+ x,. By considering 
the other quadrants we can find equations for the unit ball in each case, and 
then see that the unit ball B,, far from being an ordinary circle, must be a 
diamond: 
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What is the circumference of B,;? To be consistent we should measure B, 
with respect to our new distance measure || - ||,. Calculating the length of B, 
is particularly simple as it is composed entirely of 4 straight lines connecting 
the 4 points where the boundary of the ball intersects the axes. The length of 
the segment in the upper right quadrant is 


|1,0) — (0, I) =}, -D I} = Wy + | — 1 = 2. 


This is the same for all 4 lines so the circumference of B; is 4x2 = 8. Plugging 
this into our definition of m in (6.8), noting that as the boundary of the unit 
ball b; as “radius” 1, we have 


1 1 
Ti = [circumference of unit ball] = 78) = 4, 


The value of x implicitly depends on how we measure distance! 


2. BEYOND THE MANHATTAN NORM 


The two norms we have been considering so far — the usual Euclidean norm 
and the Manhattan norm — are just two special cases of a general family of 
norms called the p-norms. The p-norm is defined for every p € [1, 00) by 


1 
lxllp = (lal? + [xal?)"”. 


If we plug in p = 1 and p = 2 we recover equations (6.10) and (6.9) respec- 
tively. As before, this way of measuring distance has the usual nice properties. 
The triangle inequality 


IIx + yllp < lly + lly llp 


follows from the Hélder inequality, which says that if p and q are positive 
numbers satisfying 


ie 4 
—~+-=1 (6.11) 


then for any n € Nand x,y € R", 


n n Vpn Va 
3 Ixivil S (> nt?) (> pu] : 
i=1 i=1 i=1 


There’s no need to get hung up on the specifics here, just to know that || - ||, 
behaves pretty much as you'd expect distance to behave. The relation (6.11) 
says that p and g are Holder conjugate. It will reappear later! Finally, if we fix 
some point x and take the limit as p — oo we are left with the supremum 
norm, 

Il lloo = max {]x1], [x2 }- 


The unit balls of the p-norm for a selection of values of p are drawn below. 
The smallest unit ball is for the p = 1 norm, which we have seen above. The 
unit balls get monotonically larger with p until we are left with the p = oo 
ball, which is a square. 
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We wish to find the z values for these norms, but calculating the length 
of the circumference of these balls is a far more difficult challenge then pre- 
viously. We can clearly see that unless p = | or p = ov the balls aren’t com- 
posed of easy to measure straight lines. This problem prompts the more gen- 
eral question: given a way of measuring distance in the form of a norm ||- || 
and a continuous curve in the plane, how do we measure the length of that 
curve with respect to the norm? 

First, let us write the curve as a continuous function @ : [0,1] > R?. (Ev- 
erything that follows actually works for curves in R” but we will set n = 2 for 
simplicity.) We can approximate the curve a by fixing a partition 


O0O=to <t) <...<tm = 1 of [0,1], 
giving rise to a sequence of points 
Xo = Q(to),...,Xm = A(tm) 
on the curve a which we can join by a polygonal path in the following fashion: 


X1 X5 


XE = a(t) 


Our approximation of the length of the curve is the length of the polygonal 
path: 


Y> la(ti) — (a) I). 


i=1 
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To get the actual length of the curve we then take the supremum (informally 
the maximum) over all possible partitions of [0, 1]: 


length(@) = sup a(t) — w(ti-1)|| + {t:}7Lo a partition of [0, 1] 


i=1 


It is geometrically clear that this does give the correct length, but what a hor- 
rendous thing to calculate! How do even start in applying this definition of 
any of the unit balls above? 

This problem is actually very tractable. First, the continuity of a means we 
don’t need to take the supremum over all partitions of [0, 1], but it is sufficient 
to take it over equally spaced partitions. Then 


’ B34 
«(=)-a(! )| ment. (6.12) 
m m 
This is a significant simplification. 


Now let us assume that a is differentiable, with derivative a’ : [0, 1] > R?, 
and fix m € N. We can make an appeal to the mean value theorem to say that 
for every integer i € [1, m] there exists 


an Con 
Sm © a ii 


m 


length(a) = sup | Si 


i=1 


such that 


cya A) a lo()-e(Z)] 0 


This isn’t strictly speaking correct as we really have to apply the mean value 
theorem separately to each component function of a, but the general argument 
isn’t significantly altered. Now we can plug (6.13) into (6.12), giving 


m 


length(w) = sup »» ~ ||a’(s;,) || : me N 


i=1 


’ 


The quantity 


J2(sn) | 


S| 


m 
i=1 


is just a Riemann sum for the integral of the function t + ||@’(t)|| which exists 
as the function is continuous. Thus taking the supremum yields 


1 
length(a) = i} |o'(2) | ae. 
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This is a much simpler expression! Readers who have studied the differential 
geometry of curves and surfaces may recognize it; in that subject it is called 
the arc-length formula. It is interesting to see that the formula does not rely on 
the Euclidean distance function used in classical differential geometry but is a 
general formula that holds in all real normed vector spaces of finite dimension. 


3. JU IN THE P-NORMS 


We are now ready to calculate 2), the value of 7 when we use the p-norm. To 
do this we first need to parameterize B,, the boundary of the unit ball of the 
p-norm; that is, find a continuous function a : [0, 1] > R? whose image is Bp. 
We can take a shortcut by using the multiple symmetries of the unit ball; this 
way we only need to parameterize one-eighth of B,. From the diagram, 


X2 


we see that 


a(t) = (0, (i- es te [o.2-77| 


does the trick. (You can check that ||@(7)||, = 1 and hence that a(t) € Bp for 
every t.) We have 


1 
a(t) = (1.-2 Gary pt?) = (1,- (1-2)? Pt) 
P 


’ tP pl we 
wo, [1+ (Fs) | 


and hence 
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Then 
1. : 
Kp = 5 [circumference of unit ball] 


= ; [8 x length of a] 


2-1/p 1? p-l 1/p 
=4 1 —— dt. 
I Hp) 


The integral substitution u = 2t? simplifies the integral somewhat, giving 
2 1 
Se =f [wl 4 uy]? au, (6.14) 
0 


This expression is complicated, to say the least! Rather than trying to inves- 
tigate it analytically, the easiest thing to do is plot it. This involves numeri- 
cally integrating the expression (6.14) in your favourite mathematical software 
package for a selection of values of p. The result looks like this: 


> Pp 


I — ae ae “28 = = em OP - 
There is some notable behaviour. First, 2, appears to be restricted to the in- 
terval [3,4]. It has a maximum at p = 1 where 7; = 4 (this we’ve already cal- 
culated) and seems to decrease monotonically to p = 2 where it has a global 
minimum of 

Mo = mw = 3.14159... 


After p = 2 it increases monotonically. Calculating 7. in the same way we 
calculated z (recall the unit ball for 7, is an easily measured square) we find 
Too = 4. It would appear that the value of x, asymptotically tends to this 
value as p —> oo. 

Perhaps the most interesting behavior relates to the values of p for which 
Ip is the same. In the diagram we that p = 4and p = 4/3 provide an example. 
From 


we see that these p values are Holder conjugate of each other, following the 
definition in (6.11). This turns out to be a general fact: if p 4 q then we have 
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Tp = Tq if and only if p and q are Holder conjugate! Observing that 2 is the 
only number Hélder conjugate to itself, 


2 + a =1 => p=2, 
DP P 
we see that p = 2 is the only value for which zp is unique. 

One way to prove this Hélder conjugate fact would be manipulate the rep- 
resentation of zp in (6.14). As might be expected, such an endeavor is difficult 
because of the complexity of the integral representation. However, this ap- 
proach is unnecessary as it turns out that this fact is actually just a special case 
of a deeper equivalence of z values. To develop such deeper results we must 
stop restricting ourselves to the p-norms and ask: for any 2-dimensional real 
normed vector space, what can we say about its 7 value? 


4. JU IN GENERAL 


So far we have focused on generalized x for given specific distance functions, 
namely the p-norms. As modern mathematicians we should like to study the 
a concept in the abstract, and see if we determine general characteristics of it 
independent of specific norms we look at. 

In the previous section we saw that 2, was bounded between 3.14159... 
and 4. Is it bounded generally? The answer turns out to be yes: if we want the 
very intuitive triangle inequality to hold we must have 


3<n<4. (6.15) 


Why? The answer is geometric. The triangle inequality forces the unit ball to 
be convex: if the interior of the unit ball contains points x and y it must contain 
every point on the straight line joining them. This means in particular that the 
boundary of the unit ball can’t be made long by wiggling it about; for instance, 
it can’t do this: 


ee not convex 


Convexity sets an upper bound on how long the boundary can be. 

The existence of a lower bound is easier to justify. As the unit ball has radius 
of 1, any two points on the boundary of the unit ball exactly opposite each 
other will be a distance 2 apart. Going between these two points by travelling 
along the boundary of the ball will take at least as long as travelling along the 
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straight line joining them, and hence the total circumference of the unit ball 
will be at least 2 x 2 = 4. The z value must be at least 4/2 = 2. 

We have already seen that the upper bound in (6.15) is attained: in the very 
first section we found a norm, the Manhattan norm, for which z = 4. Is the 
lower bound attained? The answer is yes, and proving this is even easier, once 
you know how to draw the ball. 


A 


I 46 


by translation invariance 


L(qa;) = L(bj) = 1 


l 6 
n= 5) L(@i)=3 


i=1 


a3 ¥ 


The last general fact we will present here demonstrates the deeper reason 
that m,» = mq if p and g are Holder conjugate. If we take any normed vector 
space X — informally R” together with some norm — we can form a dual space 
consisting of the set of all linear maps from X to R. Linear means that if x, y € 
X anda, B € R then 


plax + By) = ab(x) + BP(y). 


The set of such maps, the dual space, is again a normed vector space. It has 
the same dimension of the original space, and comes automatically equipped 
with a norm induced from the original norm. In our context, this means that 
if we take any norm on the plane R? the dual space will be R? again, but with 
a different norm. We have the following fact: 


the z value of a space is exactly the same as the z value of its dual space. 


This is related to the previous section by the following theorem: the dual space 
of R? with the p-norm is precisely R? with the g-norm where g is the Hélder 
conjugate of p. This explains why numbers which are Holder conjugates of 
each other have the same zp value. 


5. IS GENERALIZED 77 GOOD FOR ANYTHING? 


Throughout this article we have seen that the value of depends on our intu- 
itive idea of distance, and hence that by looking at different ways of measur- 
ing distance we get a different value of x. In the article introduction I claimed 
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that this playing around with m has some actual mathematical application and 
wasn’t merely poking holes in one of mathematics’ most prized concepts for 
nothing. In fact, I began studying the generalized concept not for its own 
value but out of a desire to find a more geometric way of solving a particular 
kind of mathematical problem: classifying normed vector spaces. 

It is hard to describe this briefly in an accessible way, so I will just give a 
basic sketch. We have already seen the idea of a norm. An isometry between 
two normed vector spaces X and Y is a continuous bijective linear map T : 
X — Y which preserves distance: for all x € X, 


IIxllx = IIT @)lly - 


If such an isometry exists the two spaces X and Y are said to be isometric and 
are considered, from the point of view of functional analysis, to be identical. 
My project started by examining R” with the p-norms, and seeing which of 
these spaces are isometric. Trying to directly prove two different spaces are 
not isometric turns out to be quite hard, even when it is geometrically obvious 
that they aren't. So why is z relevant? 

Consider By, the boundary of the unit ball of the normed vector space X. 
If we apply our isometry to By we must get some subset of By as the distance 
value of 1 is preserved for each point. This subset must be the whole of By as 
T is surjective (onto): any point in By must be the image of some point in X, 
and this point in X must have distance 1 as T is an isometry. Thus, in symbols, 


T (Bx) = By. 
Now if we take any parameterization a : [0,1] > X of By then the function 
Toa:[0,1]J>Y 


is a parameterization of By. Then 
circumference of By = - | (7 o@)'(t))|] y de 
= i || T(@’(t)) ||, det 
=f leolear 


= circumference of By. 


The circumferences are the same, so the z values must be the same. The num- 
ber z is thus an isometric invariant: if two spaces are isometric they must have 
the same z value. An easy route to proving two spaces are not isometric is 
then to show that their z values are different. In our case this gives: if p 4 q 
and p and gq are not Hélder conjugate then R? with the p-norm is not isometric 
to R? with the g-norm. 
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LARGE DEVIATION THEORY: ESTIMATING THE 
PROBABILITIES OF RARE EVENTS 


BRENDAN WILLIAMSON 
Dublin City University 


In statistics and probability we are often interested in the “long run” behaviour 
of certain objects, be they estimators of certain quantities, random walks, 
stochastic models or any other types of random mathematical objects. In the 
most elementary of these cases, the normalised random walk 


1 n 
1= 


where X1, ..., X;, are independent and identically distributed random variables, 
the tools most widely available are the weak and strong laws of large numbers, 
and the Central Limit Theorem. The former give us a value for which we can 
be reasonably certain our random walk tends to in the long run, and the latter 
provides us with an approximate distribution for our random walk when n is 
large. Both of these methods have similar limitations. The weak and strong 
laws of large numbers don’t give any characterisation of probability relating to 
the behaviour of our random walk, and the Central Limit Theorem, although 
it attempts to, has been known to be inaccurate in calculating the probability 
of rare events. For example, when X, is Bernoulli with p = 0.2, 


P[Si09 < 1/100] + 5.09 x 107°, 


yet 

PIN (1/5, (1/25)7) < 1/100] = 3 x 1077, 
from [1]. Although both quantities are small, the Central Limit Theorem ap- 
proximation is almost 100 times the size of the actual value. (Admittedly this 
example is slightly cherry picked, partially to enable the use of standard Nor- 
mal tables, but we will return to it later.) 


1. CRAMER’S THEOREM IN R 


To remedy this discrepancy Harald Cramér proved what is now referred to 
as Cramér’s Large Deviation Theorem. It has some generalisations, and has 
inspired similar theorems, but its original form is stated as follows. 


Theorem 1. (Cramér) Let {Xn}0°_, be a sequence of independent and identically dis- 
tributed random variables with cumulant generating function A : R — (—oo, oo], 
and leh Spins 7g Ns Then 


1 
— inf A*(x) < liminf — log P[S, € B] 
xeB° noo n 


1 
< lim sup — log P[S, € B] < — inf A*(x) 
n xEB 


noo 
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for any measurable set B, where B° denotes the interior of B, B its closure, and 


A*(x) = sup{@x — A(6)}. 
OER 


This may be a lot to take in initially, especially to readers unfamiliar with 
measurability, interiors and closures, or limit inferior and superiors. The im- 
portant thing to take from this theorem is that for every set B for which P[S, € 
B] makes sense, we have an idea of the behaviour of n~! log P[S, € B] for 
large values of n. For example, let {X,}°°, be Bernoulli random variables with 


p = 0.2, as before. Then 
A*(x) = sup {ox =log(l =p pe") 
OER 
It can be shown using calculus that this supremum occurs when 


9 _ *0—p) 


pu =x) 


which only make sense when x € (0, 1). In this case 


A*(x) = x log (=) +(1—x)log (=) , 


When x < 0 or x > 1 the supremum occurs by letting 6 tend to minus infinity 
or plus infinity respectively, in which case the supremum is infinity, unless 
x = 0,1, in which case A*(0) = —log(1 — p) and A*(1) = —log p. These last 
values coincide with the expression for A* above, by defining 0 log 0 = 0. So, 
in conclusion, we have 


x 1-x 
x log (=) + (1 — x) log (=) x € [0, 1] 


+00 otherwise. 


A*(x) = 


This function is almost continuous (continuous apart from at x = 0, 1), so we 
can prove that inf,e@,5) A* (x) = infxeja,5] A*(x) for all a,b € [0, 1] and hence 
that 


1 
lim —logP[S, € (a,b)] =— inf A*(x). 
jim, 5 los PlSn € G,5)] a (x) 


Further analysis of A* will reveal that it is non-negative and has a unique zero 
when x = p = E[X,]. Also, the above statement implies that, for large n, 


P [Sn € (a,b)] + f(n) exp }-n inf, At 


for some function f that increases sub-exponentially if at all, and is possibly 
dependant on (a, b). Under this construct, returning to our original Bernoulli 
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example for p = 0.2, 


1 1 1 99 99 
P E < ml ~ f (100) exp (-100 ( 100 1°8 (=) + 790 108 (5))) 


= f (100) (1.38 x 10-8) 


which is still significantly better than our Central Limit Theorem approxima- 
tion if we assume f ~ 1. It can be proven easily that f = 1 uniformly if we 
examine P[S, = 0], so this may not be a bad assumption to make. 

Other direct observations we can make are that for two sets, say [.1, .15] 
and [.25, 3], we have 


_ log P[S, € [:1, .15]] + log PIS, € [.1,.15]] 8.3786 
lim = 1 = 


uw e 155, 
n> log P[Sn € [.25,.3]] "> 1logP[S, € [.25,.3]] 7-382 


So although both events become increasingly unlikely for large n, we can see 
that, roughly speaking, S, is more likely to deviate far below the mean than 
above the mean. 

Returning to the general statement of Cramér’s Theorem, with some re- 
strictions on the distribution of X;,a number of properties of A* can be proven. 


Lemma 1. Let X a random variable and let A* be defined as in Cramér’s Theorem. 
Then A* is a non-negative convex lower semi-continuous extended real valued func- 
tion. Moreover, if the moment generating function of X is finite in a neighbourhood 
of the origin, then X = E[X] exists as a real number and A*(x) = 0 if and only if 
x = x,and the level sets 


W(a) = {x : A*(x) <a} 
are compact for all a € [0, 00). 


Table 7.1 lists the form of A* for a number of distributions. Notice how 
the Lognormal and Cauchy rate functions don’t have unique zeros or compact 
level sets, as their cumulant generating functions are not finite in a neighbour- 
hood of the origin. 


2. LARGE DEVIATION PRINCIPLES IN GENERAL 


The applications of Cramér’s Theorem and other theorems in Large Deviation 
Theory are diverse and far reaching; they can be found, for example, in statis- 
tical mechanics, thermodynamics and hypothesis testing. However they are 
also quite technical and often esoteric, so they will not be discussed here. In- 
stead we will define and discuss the core definition in Large Deviation Theory, 
the Large Deviation Principle. First however, some preliminary definitions are 
required, mostly based on some material from functional analysis and general 


topology. 
Definition 1. If V is Hausdorff topological space, then a function I : ¥ — [0, co]isa 


rate function if I is lower semi-continuous, and its level sets Wy (a) = {x : I(x) < a} 
are compact in X. 
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Distribution Rate Function 
x log ) 

Bernoulli(p) A*(x) = +(1— x) log = x € [0, 1] 
+00 otherwise 


log (fj) -x +A xe€[0, 
Poisson(A) A*(x) = x log (=) —x x € [0, 00) 
+00 otherwise 
x log (4) 
+(n—x)log (=) 


Binomial(n, p) | A*(x) = 


—nlogn x € [0,7] 
dee otherwise 
Normal(j,07) | A*(x) = 4 (4) 
0 x >0 


Lognormal(w,o) | A*(x) = 
+co otherwise 


Cauchy(y) A*(x)=+o00VxeR 


Table 7.1: Form of the rate function A* for a number of common distributions. 


Note that in the case of Cramér’s Theorem, A* is a rate function if the un- 
derlying random variable X obeys the conditions in Lemma 1. 


Definition 2. Let X be some set, and let 2* be its power set. © C 2* is a o-algebra 
of & if 

1. & is non-empty. 

2. Aj,...,An,...€ => ig An € & (Closed under countable union). 

3. Ae => AS € J (Closed under complementation). 


Here A° denotes the complement of A. Note that by using De Morgan’s 
Laws, (ii) and (iii) can be used to prove that ¥ is closed under countable inter- 
section. (1, 4’) is referred to as a measurable space. 


Definition 3. Let (4X, t) be some topological space. Then B(4’, tr), the Borel o-algebra 
on (AX, t) is the smallest o-algebra that contains all the open sets in (Xt). If B(V) © 
» for some o-algebra X', then & is said to contain the Borel o-algebra on (XT). 


o-algebras are of huge importance in probability, especially when analysing 
probabilities on a sample space that isn’t R”. When setting up a probability 
measure on a sample space we must do so with respect to some o-algebra on 
the space, often the Borel o-algebra. With these definitions in mind we are 
now ready to define the Large Deviation Principle. 
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Definition 4. Let (V, 1) be a Hausdorff topological space, let X be a o-algebra that 
contains B(X), let {Pn[-]}92, be a sequence of probability measures on (X, X’), and 
let {an}%°_, be a positive, strictly increasing sequence tending to infinity. Then the 
pair {(Pnl-],an)}02, is said to obey a Large Deviation Principle with rate function 
I: X = (0, o] if 


1 1 
— inf /(x) < liminf — log P,[B] < limsup — log P, [B] < — inf /(x) 
xeB° n>oo Ayn n>oo 4n xeB 


for every B € XY. Ifa, =n forall n > 1 then we say {Py[-]}7-, obeys a Large 
Deviation Principle with rate function I. 


Note that by this definition, Cramér’s Theorem states that {P[S, €°]}9, 
obeys a Large Deviation Principle with rate function A* if the underlying ran- 
dom variables {X,}%-., have a moment generating function that is finite in a 
neighbourhood of the origin. Also, just like in Cramér’s Theorem, unique ze- 
ros of the rate function are sufficient to prove weak laws, regardless of the 
properties of the topological space. 


3. EXAMPLES OF LARGE DEVIATION PRINCIPLES 


There are a number of different methods of proving Large Deviation Princi- 
ples, from direct proofs, many examples of which can be found in a text by 
Amir and Dembo on Large Deviation Theory [2], to general methods of proof 
compiled by Lewis and Pfister [3]. However, we will summarise by looking at 
examples of Large Deviation Principles. 


3.1. Cramérs Theorem in R4 


Cramér’s Theorem in R¢ is simply a generalisation of Cramér’s Theorem to 
independent and identically distributed random vectors. 


Theorem 2. (Cramér’s Theorem in R@) Let {Xn}°°, be a sequence of i.i.d. random 
vectors in R4 with cumulant generating function A : R4 — (0, oo] which is finite in 
a neighbourhood of the origin. Then {P[Sn € -]}92, obeys a Large Deviation Principle 
with rate function 


A*(x) = sup {(A,x) — A()}, 


AeER® 


where (A, x) is the inner product between the vectors A and x. 


3.2. Géartner-Ellis Theorem 


The Gartner-Ellis Theorem is a generalisation of Cramér’s Theorem, where the 
assumption that {X;,}°° , are i.i.d. is dropped. We will state the theorem in the 
case where X, € R¢@ for all n. Let A, be the cumulant generating function of 
X,. The iid. nature of the sequence of random vectors is replaced with the 
following assumption. 
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Assumption 1. For each 4 € R4, 
1 
A(A) = lim —A,(nd) 
non 
exists as an extended real number. Furthermore, A is finite in a neighbourhood of the 
origin. 

Notice how this assumption holds trivially when {X,}0°_, are iid. 
Definition 5. y € R¢@ is an exposed point of A*(x) = SUP, cpa t(A, x) — A(A)} if 
forsome A € R¢ andallx 4 y, 

(A,x) — A*(y) > (A, x) — A*(x). 
Definition 6. Let Da = {A |A(A) < 00} and let DS, be its interior. A is essentially 
smooth if 

1. D, FG, 

2. A is differentiable on D‘,, 


3. limnsoo|VA(An)| = 00 for all sequences {A,}°_, in DS, converging to a 
boundary point of D%. 


Then the Gartner-Ellis Theorem is stated as follows. 


Theorem 3. (Gértner-Ellis) Let Assumption 1 hold with a function A that is es- 
sentially smooth and lower semi-continuous. Then {P[Sn € -]}?2, obeys a Large 
Deviation Principle in R4@ with rate function 


A*(x) = sup {(A,x) — AQ)}. 
AERZ 
3.3. Sanov’s Theorem 


Theorem 4. (Sanov’s Theorem) Let X bea complete seperable metric space, let {X;}?° , 
be a sequence of random variables in X with distribution . Let M() be the space 
of probability measures on X, and let t be the topology on M(%) generated by a base 


consisting of all sets of the form 
i fdv-x|<6 
x 


forx €R,5>0, f € & > R bounded and continuous. Let 


eee 
Ln =- bx, . 
n nk 


be the empirical measure corresponding to the first n observations, dy, denoting the 
unit mass at Xx. Then {P[Ln € -]}92, obeys a Large Deviation Principle with rate 
function 


Ux8.f = a 


dv 
Hiv =i dv log — 
(v||H) : 8 au 
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So in the same way that Cramer’s Theorem served to estimate probabil- 
ities of sample averages, Sanov’s Theorem estimates probabilities of sample 
densities. The topology referred to above is known as the weak topology. For 
example, if {X,}°2, are from a Normal(y, 07) distribution, and v represents 
an Exp(A) distribution. Then 


d bg aarti Se 
= eo x and dv = Ae x 
V 200 
and so 
[ese) = hewx 
Holl) = f Ae A log | ep? dx 
20 i 


=. ee —Ax (x ~~ y)? 
= Ae~** (log(Aav 2m) + ——z— —Ax)dx 
0 o 
1 fee Bo Sys SE 

= log(Ao V2) + a(7 + a2 a 8 ) 
So the probability of obtaining what “looks like” an exponential distribution 
from asample of normal random variables decays exponentially with the above 
parameter. Using calculus we can show that when y = 0,4 = 2 minimises the 
above expression, and therefore out of all exponential distributions, an empir- 
ical sample from a Normal(0, 07) is most likely to converge to one with mean 


5° 
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SIMULATING A TSUNAMI IN MATLAB 


ANTHONY JAMES MCELWEE 
University College Dublin 


1. INTRODUCTION 


This article presents a simulation that depicts a tsunami wave crossing an 
ocean. It adapts a finite-difference finite-time scheme that was implemented 
using the Lax-Wendroff method in order to solve the Shallow Water Wave 
Equations. While computational methods are time consuming and often dif- 
ficult to implement, the rewards for implementing such schemes are high and 
often visually pleasing. Indeed, there is a general feeling nowadays that prop- 
erly trained modern applied mathematicians, physical scientists, and engi- 
neers should have some understanding of numerical methods. The compu- 
tational work was coded and rendered in Matlab. This software is available 
to most undergraduate students in science and engineering disciplines. Upon 
reading this article, it is hoped the reader will be able to create an .avi film clip 
that depicts the a basic simulation of a tsunami that is at the far-from-shore 
stage of propagation. Screen grabs of such a film are pictured below. 


2. PREVIOUS WORK 


This simulation builds on previous work carried out by Moler [3]. While this 
code creates a simulation pictured below, it lacks graphical quality and is in- 
capable of being directly captured in a film file. 

Moler’s use of special loops for stability and graphics initialisation are com- 
mon coding practice but such elements make the code inaccessible to under- 
graduates who want to play with the finite-difference finite-time schemes and 
the visual parameters of the simulation. Any change could prevent the simu- 
lation from running and thus serve as a disincentive to students investigating 
numerical schemes. The new simulation created a dramatically altered code 
that is geared towards including easily adaptable graphical parameters and 
the ability to be captured to an .avi file. In order to emphasis these aspects, 
the new code abandoned much of the numerical scheme. However, a partial 
derivation of the numerical scheme from the partial differential equations will 
be outlined in the next section. 


Mathematical Scheme 


The Shallow Water Waves imitate waves whose depth is much less than the 
wavelength of the disturbance of the free surface. The model assumes that 
the vertical acceleration of the fluid is much smaller than the gravitational ac- 
celeration and that the horizontal component of the fluid velocity is uniform 
along any vertical section through the fluid (Billingham and King [1]). The 
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Figure 8.1: Screen grabs of the final simulation. 
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Figure 8.2: Screen grabs of Moler’s simulation. 
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partial differential equations used by Moler [3] correspond to those derived 
by Billingham and King [1] from the Navier-Stokes equations. These derived 
equations are 


oh rf d(uh) di d(vh) a 


ot ox oy : 
d(uh) n a(u7h + 4 gh?) 1. (uh) _ " 
ot ax oy 
d(vh) d(uvh) —0(v?h + 4 gh? 
OP OUDI) 4, Bele) 2h 
ot ox Ox 


where h corresponds to the height of the water and g is the gravitational ac- 
celeration. The two-dimensional velocity field is denoted in terms of u and 
v and their products with h are proportional to the water’s momentum. The 
perspective taken is that the mass and momentum of the water is conserved 
at all times. A number of numerical difference methods could be considered 
that are suitable for conservation problems but the Lax-Wendroff method is 
satisfactory for this situation. 

The next step is to recast the equations so that they are in the form of a 
hyperbolic partial differential equation which obeys the conservation laws. To 
do this, the following three vectors were formed: 


h uh vh 
U=|uh|, FU)=|u2h+ sgh? , GUW)= uvh , 
vh uvh v*h + $gh? 


which gives our model equation 


aU, OF) | AGU) _ 


0. 
ot ox oy 


The U corresponds to the height matrix, while Fis the matrix that stores the 
zonal velocity and G stores the meridional velocity. 


Numerical Scheme 


Unfortunately, any land masses or irregular geometry have to be discounted 
in order to keep the scheme relativity simple. This implementation is over a 
square grid with regular divisions. Such a grid is the easiest to implement in 
Matlab. The grid is also stepped so that time and space are subdivided, as 
pictured in Moler’s brief guide to the original code. 

Such a grid can be revealed by altering the code provided later in article so 
that Edge_Alpha has a value close to 1, rather than 0. 

As already stated, the original code uses the Lax-Wendroff numerical 
scheme to solve the new model hyperbolic partial differential equation. The 
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Figure 8.3: The grid showing the values for the vectors positioned in a half 
step position on both axis. 


steps involved in the scheme are described in most computational fluid dy- 
namics books and freely available on the internet. To quickly summarise, the 
finite difference representation of the model equation is derived from the two 
dimensional Taylor series expansion of the dependant variables. The time 
derivative of the model equation is taken to obtain a new equation and then 
the model equation and the new equation are substituted back into the Taylor 
expansion of the dependant variables. In one direction, say x to illustrate, the 
equation becomes 


uttl un at Ui. — iy 4 LAL? Uj, — 2U?' + Uj", 
: : 2Ax 2 (Ax)? 


Moler goes further and half steps the time and space steps. 


The full set of equations for the two time steps are available in Moler [3]; 
the present author sees no benefit in rehashing them here. At this point we 
depart from Moler and decide to take action with a view to improving the 
visual quality of the simulation. For the simulation will assume that only the 
y-components of the equations matter and set the x-components equal to the 
y-components. This shortens the code and since this isn’t a real tsunami no 
one should really care! The code could be rectified by the reader to fully reflect 
the nature of the Shallow Water Waves. Provided the stability of the numerical 
scheme is maintained and the film looks like a wave propagating across a large 
body of water, the goal of this exercise will have been achieved. This allows 
us to get into the visual aspect of the simulation. 
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3. METHODOLOGY 


The adapted code is now explained in the order that it appears in the pro- 
gramme. The reader should try adapting the code parameters further and 
experiment with the visual effects. If any of the code remains unclear after 
reading the explanations, extra help can be found using the help function in 
Matlab. 

The first thing to do is to declare the function in a new .m file that is saved 
in the Matlab workspace. This allows the code to be easily updated. Here 
the function is called FDFTSWWIUMM, (Finite Difference Finite Time Shallow Wa- 
ter Waves Irish Undergraduate Mathematical Magazine). Clearing commands 
are included so that when the simulation is restarted, the old values will not 
interfere with the fresh simulation. 


function FDFTSWWIUMM 


Ht 

2 

3 | 3% Clear MATLAB 

4 |clear all;clc;clf 


The next section establishes some basic scene settings. The axes are stripped 
away and the simulation freezes its aspect ratio, otherwise the simulation will 
lack any sort of aesthetic appeal. The z-axis, corresponding to the height of 
the water, is also fixed. A failure to set the z-axis would result in the scene 
having a jittery appearance, even if the amplitude of the perturbation wave is 
relatively small. 

Setting the handle for the background will probably be new to the reader. 
This gcf is setting the background colour using the RGB colour code with a 
base of 255 integer values. The surface handle is set in a similar manner. The 
water_range generates a vector for the surface colour. 


%% Scene Settings 

axis vis3d off 

zlim([-50 50]) 

set(gcf,’ color’ ,[57/255 176/255 121/255]); 
water_range = [20/255 80/255 118/255]; 


Oo OND MH 


The next step is to set any constants such as gravity, the time step and the 
number of time steps the simulation will last. Initially, the time length is set at 
40. It is kept short since a requirement to change some basic parameter can be 
spotted by the reader in that time frame, removing the need to repeatedly run 
ling simulations. 


10 |%% Constants 

11 | grav = 9.81; 

12 |dt = 0.5; 

13 | sim_steps=40; 


The next section declares a batch of variables that are the handle property 
settings for surfaces. This is not a complete list of possible variables, but they 
are comprehensive for this type of simulation. They will define the handle’s 
set properties. In order to make it easier to adapt the variables, they have 
been named in a similiar fashion to those property titles. This allows you to 
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leave the surface handle commands alone and therefore reduce the chances of 
accidentally introducing an error into the code. 


14 | %% Set Handle Property Settings. 
15 | Ambient_Strength=0.5; 

16 | BackFace_Lighting=’1lit’; 

17 | Diffuse_Strength=0.8; 

18 | Edge_Alpha=.1; 

19 | Edge_Color=water_range; 

20 | Edge_Lighting=’ phong’ ; 

21 | Erase_Mode=’normal’ ; 

22 | Face_Alpha=0.9; 

23 | Face_Color=water_range; 

24 | Face_Lighting=’ phong’ ; 

25 |Mesh_Style=’both’; 

26 | Normal_Mode=’auto’; 

27 | SpecularColor_Reflectance=.5; 
28 | Specular_Exponent=10; 

29 | Specular_Strength=.5; 


The section above offers a lot of possibilities and is the key to enhancing the 
visual properties of the simulation. The obvious place to start is by varying the 
numbers but some of the settings are set by string commands. More informa- 
tion about the other available options for these string commands are available 
at www.mathworks.co.uk/help/matlab/ref/surface_props.html. For exam- 
ple, gouraud shading could be used instead of the phong settings. These com- 
mands have a vast amount of physical attributes that go to the core of com- 
puter graphics. By varying these properties, the simulation can take on a new 
appearance. 

The next section is a straightforward creation of the square, uniformly di- 
vided grid. Instability will occur if the grid step and the time step are beyond 
certain values. This is manifested in a simulation that fails to run, or blows 
up rapidly. Usually, the parameters can be choosen before running the simu- 
lation. For example, the Lax-Wendroff scheme usually has the stability condi- 
tion that corresponds to the time step over the spatial step having a ratio less 
than one. Since a lot of the numerical scheme has been dumped, no definition 
of the stability for the new situation exists. The reader should experiment by 
changing these variables to understand the significance of selecting step sizes. 


30 | 3% Grid Parameters 


31 |X_length = 100; 

32 |Y_length = X_length; 
33 |X_n = 35; 

34 |Y_n = X_n 

35 |dx = X_length/(X_n-1); 
36 | dy = dx; 


37 
38 | 3 Grid Generation 

39 | [x_grid] = 0O:dx:X_length; 
40 | [y_grid] = O:dy:Y_length; 


The next section sets up the initial conditions of the surface so that the 
ocean is calm and flat. One could add some noise to the surface to mimic 
the normal action of the sea but due to the scale that tsunamis are modelled 
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on, such noise based disturbance would be pointless. There were boundary 
conditions included in the original code that were perfectly reflective and so 
that the waves bounced off the sides of the grid. This is not the behaviour of 
a real tsunami and since some of the numerical scheme was thrown away the 
author has decided to get rid of the boundary conditions too. The reflective 
nature of the boundaries in the new code is attributed to the butchered scheme 
outlined later in the article. 


41 | %% Initial Conditions 
42 |height = zeros(X_n,Y_n);height_x = zeros(X_n,Y_n);height_y 
= zeros(X_n,Y_n); 


43 | zonal = zeros(X_n,Y_n);zonal_x = zeros(X_n,Y_n);zonal_y = 
zeros (X_n,Y_n); 

44 |meridional = zeros(X_n,Y_n);meridional_x = zeros(X_n,Y_n); 
meridional_y = zeros(X_n,Y_n); 


Now the simulation needs a perturbation on the surface so that we can 
get some action into the simulation. A tsunami-like disturbance can be cre- 
ated which propagates in time and space by plotting a simple ridge uniformly 
across the width of the surface. The author recommends that the reader tries 
different plots and perturbations. For example, the use of a circle could mimic 
the impact of an asteroid hitting the ocean. 


45 | %% Pertubation Of Surface [Tsunami Wave] 

46 | for i=1:X_n 

47 Ts = 1; 

48 wavelength = 1000; 

49 amplitude=5; 

50 phase=((((X_n-1)/(wavelength) ) -(dt/Ts))+((CY_n-1)/( 
wavelength) )-(dt/Ts))); 

51 height (i,15)=amplitude*sin (2*pi*(1/Ts)+phase) ; 

52 | end 


Before the time progression of the simulation is implemented, Matlab needs 
to be able to capture and export the simulation in an .avi file. The full path is 
set to where the simulation will be exported. The gca command shown is 
included so that the axis are returned to each frame captured by Matlab. If 
the gca command is omitted, the axis settings will default every time the loop 
executes and Matlab will stall almost immediately. 


53 | 3% Movie Capture Part0Ol 

54 |vidObj = VideoWriter(’C:\Desktop\latestil.avi’); 
55 | open(vidObj); 

56 

57 | set(gca,’nextplot’,’replacechildren’) ; 


The numerical part of the programme, which will dictates how the surface 
behaves on the screen, will now be dealt with briefly. The scheme is looped 
for the designated number of time steps. For each time step, all of the matrices 
are recalculated using the previous values and the following scheme: 
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2dt* 
Ay(,j) = 2HG¢,j) — Hx@,j) + Gedy {Hi+1,/) 


+AG-1,) + Hoj+y + Ha si-y — 4HG,p} > 
dt 
Fy@,j) = Fxa,j) - 87 (A417) — AG-1,;)), 
dt 
GyG,j) = GxG,j) — 8 ay Heit — Aq,j-1)- 


The loop sets the next step values of the x-axis values to the current height 
values and then sets the next step values of the height values to the current 
y-axis values. The result holds no real mathematical significance. It is merely 
a happy accident but the scheme works for our purposes. 


58 | 33 Time Progression 
59 | for n=1:sim_steps 


60 for i=2:X_n-1 

61 for j=2:Y_n-1 

62 zonal_y(i,j)=zonal_x(i,j)-grav*(dt/dx)*(height (i+1,j) 
-height(i-1,j)); 

63 meridional_y(i,j)=meridional_x(i,j)-grav*(dt/dy) *( 
height (i,j+1)-height(i,j-1)); 

64 height_y (i,j) =(2*height (i,j)-height_x(i,j)+(((2*dt72/ 


dx*dy))*(height (i+1,j)+height (i-1,j)+height (i, j 
+1) +height (i,j-1) -4*height(i,j)))); 

65 end 

66 end 

67 height_x=height ; 

68 height=height_y; 


The next section gives the code for the surface handle and its properties. 
Unless the reader wants to insert a new property, it is best to leave this section 
alone and just change the variables previously declared. 


69 |surf_handel = surf(x_grid,y_grid,height) ; 

70 set(surf_handel,... 

71 >AmbientStrength’,Ambient_Strength,... 

72 >BackFaceLighting’,BackFace_Lighting,... 

73 ?DiffuseStrength’ ,Diffuse_Strength,... 

74 >EdgeAlpha’,Edge_Alpha,... 

75 >EdgeColor’,Edge_Color,... 

76 >EdgeLighting’,Edge_Lighting,... 

77 »EraseMode’,Erase_Mode,... 

78 ?>FaceAlpha’,Face_Alpha,... 

79 >FaceColor’,Face_Color,... 

80 >FaceLighting’,Face_Lighting,... 

81 >MeshStyle’,Mesh_Style,... 

82 ?>NormalMode’ ,Normal_Mode,... 

83 >SpecularColorReflectance’, 
SpecularColor_Reflectance,... 

84 >SpecularExponent’,Specular_Exponent,... 

85 >SpecularStrength’ ,Specular_Strength) 


This final part defines and positions the global lighting source for the scene. 
The generated frame gets pushed to the video object where it is stored and 
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finally closed. Now run the whole programme from your script environment 
and it’s lights, camera, action! 


86 | light(’Position’,[100 40 80],’Style’,’infinite’); 
87 
88 3% Camera Position 
89 campos (’auto’) 

90 
91 | 3% Movie Capture Part02 


92 currFrame = getframe; 
93 writeVideo(vid0bj ,currFrame) ; 
94 end 


95 
96 | 3% Movie Capture Part03 
97 | close (vidObj) ; 


4. CONCLUSIONS 


A simulation has been produced that looks like a tsunami wave propagating 
across an ocean, even if the readers know that it is not a proper representation 
of the Shallow Water Waves. By varying the properties, the reader can alter the 
visual output of the scene and change how the scene is rendered by Matlab. 
Some simple examples are shown below. 


The author recommends Moler [3] for some further information on the 
original code that provided the ground work of this article. A future project 
could be to tackle a full implimentation of the Shallow Water Waves using Lax- 
Wendroff, or using the MacCormack method, in combination with the surface 
handle and .avi export options. The use of a moving camera should also be a 
investigated. Those interested in more advanced material in the area of mod- 
elling large sea surfaces should consult the freely available Loset [2]. 
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Figure 8.4: (a) Increased initial amplitude. (b) Reducing Face_Alpha and in- 
creasing the Ambient Strength. (c) Removing the cell faces and plotting the 
mesh. 


COMPUTATIONAL ANALYSIS OF THE DYNAMICS OF 
THE DIPPY BIRD 


GLENN MOYNIHAN AND SEAN MURRAY 
Trinity College Dublin 


A simple, quantitative, model was constructed and solved in Mathematica that 
replicated the short-, medium- and long-term motion of the Dippy Bird (DB). 
An extended model for two coupled DBs was then constructed and produced 
motion that displayed stable in-phase and anti-phase motion as well as pe- 
riod doubling for the second bird. Bifurcation diagrams were produced in an 
attempt to observe chaotic behaviour for which we found some signs. 


1. INTRODUCTION 


The DB is an entertaining toy that has adorned the desks of enthusiasts since 
its invention in the early 20 century. Much of its appeal comes from its appar- 
ent display of perpetual motion; however, upon a closer look we can see that 
the bird is in fact a thermal engine that is obeying simple thermodynamical 
principles. 

The source of the DB’s power comes from the temperature difference in- 
duced by the evaporation of water on the bird’s head. Figure 9.1 is a schematic 
diagram of the bird outlining the basic mechanism. The bird consists of two 
glass bulbs (head and body), connected by a thin glass tube (neck) that is sub- 
merged in a volatile liquid in the body. The liquid and its vapor are all that 
reside within the bird - there is no air present. The glass ensemble is mounted 
on a pivot that allows it to rotate backwards and forwards and the entire sys- 
tem is in thermal equilibrium with the environment. 

The bird begins tilted at an angle (near the upright position) and when 
released starts to oscillate. The head is covered in felt that retains water. If 
the humidity of the environment is less than 100%, spontaneous evaporation 
of the water occurs. As the water evaporates it absorbs the latent heat of va- 
porisation from the head, thereby lowering the temperature of the head. This 
decrease in temperature causes the pressure of the vapor in the head to drop 
according the Clausius-Clapyron relation. The pressure drop forces the liquid 
to rise up the neck thus raising the bird’s centre of mass. As the centre of mass 
transitions to the head the amplitude of the oscillation decreases and the bird 
simultaneously begins to fall further forward with each swing. At a critical 
angle (near the horizontal position) an air channel forms in the neck allowing 
the pressure between the head and body to equalise (image on the right in fig- 
ure 9.1). The liquid rushes back to the body, lowering the centre of mass and 
beginning the next cycle. 

The DB possesses three intrinsic time-scales: the first (short) time-scale is 
the time taken for each single oscillation, the second (medium) timescale is that 
between successive “dips” where the cycle resets and the third (long) time- 
scale is one that demonstrates the increase in the cycle period. The rate of 
evaporation of the water is proportional to the mass of water present on the 
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UU 


Figure 9.1: Schematic diagram of the DB in the upright and horizontal posi- 
tions. In the right-hand image we note the protrusion of the neck from the 
liquid allowing for pressure equalisation. 


head. As the water diminishes its evaporation rate slows so the whole cycle 
takes longer to complete. We aim to reproduce each of these time-scales in our 
model. 


2. MOTIVATION 


The motivation behind this project was to arrive at a simple, qualitative model 
of the DB that would describe the general physics of motion without the need 
for an exhaustive approach. Our initial results will be compared to that of two 
papers (Lorenz [1] and Guemez [2]) that investigated the DB more extensively 
in order to highlight the successes and shortcomings of our own model. 

Lorenz [2] constructed a detailed computational model based on the ther- 
modynamics of the bird and acquired the relevant physical parameters from 
experiment. Guemez [1] employed an experimental approach and using elec- 
tronic devices measured the dependence of the period of the DB on environ- 
mental conditions. 


3. OUR MODEL 


Single Dippy Bird 


The DB is, in essence, a damped, compound pendulum whose mass varies 
with time due to the migration of an internal liquid. However, in order to 
preserve the simplicity of the project we omitted the complicated evolution of 
the internal liquid in favour of a time-dependent, external driving-force that 
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would substitute nicely for the influence of the changing mass, torque and 
moment of inertia that would have made our model grossly more complicated 
and similar to that of Guemez [2]. 

We constructed a model based on a damped, driven pendulum and, by 
equating the torques arising from rate of change of angular momentum, the 
mass of the bird, the damping force and the driving force, we arrived at our 
equation of motion 


d? d 

We (m1? 6(t)) + mgl sin(6(t)) + bl’ (0) —IF(t,0) =0 (9.16) 

where m is the mass of the pendulum, / the (constant) string length, b the 

damping coefficient, 9 the angle between the string and the vertical and F is 
the external driving force. 

We now move to non-dimensionalise the equation by introducing a dimen- 


sionless time variable 
T=T,/ = wol. 


fa 4ar_ fee d> _ g d? 
dt dtdt Vidt’ dt2 1 dt’ 


Subbing these into equation (9.16) we arrive at 


Then 


2 F(t,6 
Oe eyo 22 NEN) 2G 


9.17 
dt? Mo at mg ow) 
We now define two control variables 
F 
pw FEB) bas, Des, 
mg mag 
yielding 
6+ sin(9) + BO— f(t, 9) = 0. (9.18) 


where the driving force, f, was tailored such that it would replicate the influ- 
ence of the elements we chose to neglect. It was determined a posteriori and set 
to be 


0, 0 ca Omax . 


f(t.0) = dete eee (9.19) 


It takes the form of a sawtooth function that increases linearly until the maxi- 
mum angle is reached where it resets to zero. 

Including an increasing force ensures that, initially, the pendulum is al- 
lowed oscillate as normal but as time progresses the mass is pushed towards 
higher angles. Once it reaches the maximum angle the force suddenly cuts 
off (simulating the equalisation of pressure in the DB and the lowering of the 
centre of mass) thus restarting the cycle. 
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As time progresses the rate of increase of force, a, should diminish. This is 
to simulate mass evaporation in the physical DB and, in order to comply with 
this effect, we set a to decay exponentially with time; thus 

ave Pt, 
where p < | is our third input parameter. 

We tested our model for various input parameters. Equation (9.18) was 
solved using the NDSolve function in Mathematica and values of 6, w and f 
were recorded at each step in order to produce plots of the motion. 


Coupled Dippy Birds 


We extended our previous model to that of two coupled, damped, driven har- 
monic oscillators described in figure 9.2. 


Figure 9.2: Sketch diagram of coupled DB system with the relevant physical 
parameters labeled. The origin is taken to be the position of Bird 1 mass in 
equilibrium. 


We derived the following equations of motion: 
ue ‘ 
qa liu) + mygly sin(61) + dO, = 1, Fy(t, 1) = cl, Ad cos(6; = od) => 0, 


a : 
aa (2) + m2zgl, sin(02) + b262 — 12 F2(t, 02) + cla Ad cos(¢ — 92) = 0, 


where all the symbols have their usual meanings: c is the coupling constant, 
Ad is the displacement of the spring from equilibrium and ¢ is the angle of 
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the spring with respect to the horizontal. The last term in each equation is 
the torque produced by the spring. The rest length of the spring is defined 
as the distance between the two masses as they hang in stable equilibrium, 
do = VX2 4 Y2 . 

Non-dimensionalising our new equations in a similar manner to before we 
arrive at 


6; + sin(O,) + B16; — a1 — y1 Ad cos(@ — ¢) = 0, (9.20) 
> + sin(42) + B22 — a2 — y2Ad cos(# — 62) = 0, (9.21) 
where P 
yi= fori = 1,2. 
&mMj 


Equations (9.20) and (9.21) were simultaneously solved using NDSolve and cal- 
culations of the same parameters were made as in the single DB case. 


4. RESULTS: SINGLE DIPPY BIRD CASE 


We now demonstrate the results of our model by comparing it to those ob- 
tained by Lorenz [1] and Guemez[2]. 


Angular Evolution 


Our first goal is to highlight the qualitative aspects of our plots in figures 9.3a 
and 9.3b that compare to the computational-based figure 9.4, taken from [2], 
and the experimental-based figure 9.5 taken from [1]. 


o We first notice the effect of damping leading to a decrease in amplitude 
with each successive oscillation of the DB. This comes as no surprise as, 
for the first few oscillations, the force is quite small. This effect is also 
seen in figure 9.4. 


fe} 


We next observe the shift in the equilibrium position as the force be- 
comes dominant. This drives the angle of the bird to higher angles (or 
in the case of Lorenz, lower angles). 


fe} 


The third, less apparent, trait is the lengthening of the period of oscilla- 
tions; i.e., the third time scale. To highlight this effect we set p = 0.002 
and plotted three successive cycles in Figure 9.3b. When compared to 
Figure 9.5b we see the general structure of the plots is quite similar, how- 
ever our period lengthening is much more obvious. 


Increasing Period 


We then plotted the period of oscillations versus the cycle number for 14 suc- 
cessive cycles. The results are plotted in figure 9.6. We see that the slope is 
quite linear in contrast to figure 9.7 of Guemez [2] where we see relatively 
flat curve and then a sharp rise. We feel the general behaviour is comparable, 
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(b) 


Figure 9.3: Angular evolution of DB. In both plots the x-axis marks time and 
the y-axis angular displacement of the bird from the vertical position of rest. 
Figure (a) has p = 0 and shows two stable cycles while figure (b) has p = 
0.002 and shows three cycles with increasing period. We thus note the high 
sensitivity of the model on the parameter p. 


6 /rad 


Figure 9.4: Guemez et al. Computational simulation of DB angle vs time. 
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Figure 9.5: Lorenz experimental measurement of one cycle in (a) and three 
cycles in (b) of the DB using electronic circuitry, with the angle taken with 


respect to the horizontal. 
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Figure 9.6: Our model: dip period vs time showing a strong linear relationship 
for p = 0.002; to be compared with figure 9.7. 
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Figure 9.7: Guemez: dip period vs time measured in an enclosed chamber 
using a digital balance. 


however, since for p < 1 we can approximate: e?4* ~ 1 + pAr and arrive at 
a linear relationship. 

The equation of the line given by data analysis is y = 36(1 + 0.002x) which 
fits very well to the linear approximation given. We would therefore assume 
that, over a longer time, we would have acquired the results in figure 9.7. 


Phase Portrait 


From the phase space diagram plotted in Figure 9.8 we get the most complete 
picture of the physics of our system. We can immediately see that it is a self- 
repeating cycle given that each cycle traces the same path in phase-space. We 
can also see that it is a dissipative system: the radius of the path decreases with 
each rotation while we can also clearly see the shift in equilibrium towards the 
later end of the cycle. 

We can then conclude that our system should not conserve energy. Indeed, 
the physical DB does work and thus does not conserve energy itself. 


Further Analysis 


We now discuss the dissimilarities. The most noticeable difference is that of 
the gradient of the envelope enclosing our oscillations. This contrast, however, 
is quite minute and does little to diminish the quality of our model with respect 
to the actual physical process. This can obviously be attributed to the over- 
simplification of our model and is, arguably, a trivial disparity to encounter in 
light of the numerous successes. 

We also note that our increasing period plot in figure 9.6 did not match 
that of figure 9.7. However this is not of great importance as figure 9.7 was 
measured in an isolated, controlled environment under specific conditions; we 
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Figure 9.8: Phase portrait of DB for three cycles sans decaying force. We clearly 
see the stability of the model as the paths are all identical. 


should therefore only concern ourselves with the general behaviour, which we 
have acquired. 


5. RESULTS: COUPLED DIPPY BIRDS 


We now end our discussion of the results of Guemez and Lorenz, having es- 
tablished a satisfactory comparison, and discuss the unique results obtained 
by this project. Having assembled our model as discussed in section 3, we 
tested our model for simple in-phase and anti-phase cases to ensure we would 
achieve expected behaviour. Once satisfied we then began to vary a in order 
to observe period doubling and/or tripling, and signs of chaos. 


In-phase and Anti-phase 


Since the focus on this part of the project was on deriving chaotic behaviour 
we set p1 = p2 = 0 for simplicity as we would not expect the force decay to 
be influential in the coupled case over short time periods. This also shortened 
our computing time. Hence both birds have identical conditions and from the 
overlayed phase-portrait in Figure 9.9a, we see the DBs’ paths in phase-space 
are identical, which is the expected result. 

We then set the DBs such that they should oscillate in anti-phase. The re- 
sultant over-layed phase-portrait in figure 9.9b clearly shows the anti-phase 
behaviour. 
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Figure 9.9: Overlayed phase portraits of birds 1 and 2. 9.9a) Birds are in phase, 
tracing the same path in phase-space and in 9.9b) in anti-phase for three cycles 
having paths out of phase by 180 degrees. 
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Period Doubling, Tripling 


Having established a satisfactory model we wished to observe period doubling 
or tripling in bird 1 by varying a2. We significantly reduced our step size to 
j = 1077 to increase our accuracy. Figures 9.10a and 9.10b show the periods 
of birds 1 and 2 respectively for each cycle. We can clearly see that after an 
initial settling period of a few cycles the birds assume a repeated pattern with 
bird 2 displaying period doubling. 


Perkxid! 


5 10 15 x 


(b) Bird2 Period 


Figure 9.10: Periods of DBs demonstrating in 9.10a), regular behaviour for bird 
1 and in 9.10b) period doubling for Bird 2, after some initial settling time. 


CHAOS 


Attempts were made at establishing bifurcation diagrams of period vs a for 
both birds and observe chaotic behaviour. We varied a, from 0.01 to 1 in in- 
crements of 0.01. The diagrams obtained are given in Figures 9.11a and 9.11b. 
It is clear these are not typical bifurcation diagrams, but figure 9.11a does dis- 
play hints of chaotic behaviour. However this is not a convincing example and 
more analysis is needed before we can say for sure. The apparent randomness 
of the period may indicate that chaotic behaviour is present but this may be a 
symptom of the aforementioned relaxing period. 
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Figure 9.11: Bifurcation Diagrams for Bird 1 in 9.11a, and Bird 2 in 9.11b) which 
show some signs of chaotic behaviour. 


6. CONCLUSIONS AND DISCUSSION 


Single Dippy Bird Model 


Our single bird model was quite successful in replicating the general motion 
of a physical dippy bird and compared very well with the results of Lorenz 
[1] and Guemez [2]. We were able to devise and observe all three time scales 
thus completing our physical description. Given the simplicity of our model, 
we believe further work into fine-tuning it based on any physical observations 
would be extraneous and unnecessary as it provides a thorough and quanti- 
tative representation as it is. 


Coupled Bird Model 


For the coupled birds we validated our model by demonstrating the predicted 
in-phase and anti-phase motion and went further to observe period doubling 
for Bird 2. We do however note that period tripling was observed for simu- 
lations where the time-step was of order 107? but then collapsed into period 
doubling where in the same simulation the time-step was of order 10~’. We 


2013] DYNAMICS OF THE DIPPY BIRD 85 


concede that the time-step in our observed period doubling may not have been 
fine enough to avoid this error. 
Hints of chaotic behaviour were observed in our simulations, however we 


Physical Model and Future Work 


Although no future work will be done on this project, had we more time, we 
would have liked to investigate the effects of a non-linear force model for the 
single DB model. We would also have liked to test our model against a real 
DB and try to make some predictions as to its motion. 

We would also have liked to re-run the simulations in the period doubling 
case with as small a time-step as possible to erase any doubt that it was indeed 
present, given that period tripling was seen to collapse when a smaller time- 
step was used. 

More simulations would certainly have been carried out on the chaotic 
model to fully exhaust any avenues in which it might be observed. We would 
take the variable parameter over a longer range, or vary $1 or y; also to see 
their influence. However, we can not say at this time if chaotic behaviour will 
be determined for certain. 
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